Tugas Praktikum Sains Data
library(tidyverse)
library(splines)
library(ISLR)
library(dplyr)
library(ggplot2)
library(SpatialExtremes)
library(caret)
library(rsample)
library(ggpubr)
library(cowplot)head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name
## 1 chevrolet chevelle malibu
## 2 buick skylark 320
## 3 plymouth satellite
## 4 amc rebel sst
## 5 ford torino
## 6 ford galaxie 500str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...Soal 1
Peubah respon yang dipilih adalah mpg dan peubah-peubah prediktor yang dipilih yaitu displacement, weight, dan horsepower.
Plot data mpg- displacement , mpg- weight, mpg- horsepower
3 peubah prediktor dengan peubah respon mpg tidak memiliki hubungan linear dikarenakan semakin besar nilai displacement, weight dan horsepower diketahui bahwa semakin kecil nilai mpg dan membentuk kurva yang menurun (memungkinkan adanya nilai mpg negatif).
Soal 2 : mpg - displacement
Regresi linear
model_linearA=lm(mpg~displacement,data=Auto)
summary(model_linearA)##
## Call:
## lm(formula = mpg ~ displacement, data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9170 -3.0243 -0.5021 2.3512 18.6128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.12064 0.49443 71.03 <2e-16 ***
## displacement -0.06005 0.00224 -26.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.635 on 390 degrees of freedom
## Multiple R-squared: 0.6482, Adjusted R-squared: 0.6473
## F-statistic: 718.7 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.50, color="red") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = F)+
theme_bw() 
set.seed(1234567890)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_linearA <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ displacement,
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_linearA## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.46 3.38
## 2 4.69 3.34
## 3 4.58 3.84
## 4 4.26 3.16
## 5 3.95 3.05
## 6 5.24 4.21
## 7 6.35 4.44
## 8 4.71 3.71
## 9 3.91 3.14
## 10 3.60 2.81mean_metric_linearA <- colMeans(metric_linearA)
mean_metric_linearA## rmse mae
## 4.575194 3.509608Regresi polinomial (visual)
- Regresi polinomial derajat 1
model_polinomialA = lm(mpg~poly(displacement,1,raw = T),data=Auto)
summary(model_polinomialA)##
## Call:
## lm(formula = mpg ~ poly(displacement, 1, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9170 -3.0243 -0.5021 2.3512 18.6128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.12064 0.49443 71.03 <2e-16 ***
## poly(displacement, 1, raw = T) -0.06005 0.00224 -26.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.635 on 390 degrees of freedom
## Multiple R-squared: 0.6482, Adjusted R-squared: 0.6473
## F-statistic: 718.7 on 1 and 390 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="red") +
stat_smooth(method = "lm",
formula = y~poly(x,1,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
- Regresi polinomial derajat 2
model_polinomialA2 = lm(mpg~poly(displacement,2,raw = T),data=Auto)
summary(model_polinomialA2)##
## Call:
## lm(formula = mpg ~ poly(displacement, 2, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2165 -2.2404 -0.2508 2.1094 20.5158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.204e+01 1.078e+00 39.011 < 2e-16 ***
## poly(displacement, 2, raw = T)1 -1.379e-01 1.113e-02 -12.388 < 2e-16 ***
## poly(displacement, 2, raw = T)2 1.685e-04 2.366e-05 7.122 5.17e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.365 on 389 degrees of freedom
## Multiple R-squared: 0.6888, Adjusted R-squared: 0.6872
## F-statistic: 430.5 on 2 and 389 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="red") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
- Regresi polinomial derajat 3
model_polinomialA3 = lm(mpg ~ poly(displacement,3,raw = T),data=Auto)
summary(model_polinomialA3)##
## Call:
## lm(formula = mpg ~ poly(displacement, 3, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.6791 -2.3900 -0.2987 2.1156 20.4528
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.445e+01 2.591e+00 17.157 < 2e-16 ***
## poly(displacement, 3, raw = T)1 -1.769e-01 3.975e-02 -4.451 1.12e-05 ***
## poly(displacement, 3, raw = T)2 3.452e-04 1.743e-04 1.980 0.0484 *
## poly(displacement, 3, raw = T)3 -2.363e-07 2.310e-07 -1.023 0.3069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.365 on 388 degrees of freedom
## Multiple R-squared: 0.6896, Adjusted R-squared: 0.6872
## F-statistic: 287.4 on 3 and 388 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="red") +
stat_smooth(method = "lm",
formula = y~poly(x,3,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
- Regresi polinomial derajat 4
model_polinomialA4 = lm(mpg~poly(displacement,4,raw = T),data=Auto)
summary(model_polinomialA4)##
## Call:
## lm(formula = mpg ~ poly(displacement, 4, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.7755 -2.3666 -0.2723 2.1005 20.4053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.536e+01 5.752e+00 7.887 3.2e-14 ***
## poly(displacement, 4, raw = T)1 -1.966e-01 1.177e-01 -1.670 0.0956 .
## poly(displacement, 4, raw = T)2 4.858e-04 8.121e-04 0.598 0.5500
## poly(displacement, 4, raw = T)3 -6.388e-07 2.282e-06 -0.280 0.7797
## poly(displacement, 4, raw = T)4 3.965e-10 2.237e-09 0.177 0.8594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.37 on 387 degrees of freedom
## Multiple R-squared: 0.6897, Adjusted R-squared: 0.6865
## F-statistic: 215 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="red") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
- Regresi polinomial derajat 5
model_polinomialA5 = lm(mpg~poly(displacement,5,raw = T),data=Auto)
summary(model_polinomialA5)##
## Call:
## lm(formula = mpg ~ poly(displacement, 5, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.3360 -2.3445 -0.2895 2.1635 20.3439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.675e+01 1.279e+01 2.874 0.00428 **
## poly(displacement, 5, raw = T)1 4.084e-02 3.360e-01 0.122 0.90331
## poly(displacement, 5, raw = T)2 -1.890e-03 3.252e-03 -0.581 0.56148
## poly(displacement, 5, raw = T)3 1.019e-05 1.454e-05 0.701 0.48363
## poly(displacement, 5, raw = T)4 -2.238e-08 3.028e-08 -0.739 0.46017
## poly(displacement, 5, raw = T)5 1.791e-11 2.373e-11 0.754 0.45103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.373 on 386 degrees of freedom
## Multiple R-squared: 0.6901, Adjusted R-squared: 0.6861
## F-statistic: 171.9 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="red") +
stat_smooth(method = "lm",
formula = y~poly(x,5,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
Regresi polinomial (empiris)
- Regresi polinomial berderajat 1
library(rsample)
set.seed(123)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyA1 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(displacement,1,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_polyA1## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.12 3.28
## 2 5.91 4.33
## 3 3.81 2.94
## 4 4.18 3.50
## 5 4.16 3.01
## 6 4.68 3.31
## 7 5.72 4.53
## 8 5.05 4.07
## 9 4.83 3.60
## 10 3.17 2.55menghitung rata-rata 10 folds
mean_metric_polyA1 <- colMeans(metric_polyA1)
mean_metric_polyA1## rmse mae
## 4.562513 3.511871- Regresi polinomial berderajat 2
library(rsample)
set.seed(1234)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyA2 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(displacement,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_polyA2## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.78 3.70
## 2 5.10 3.86
## 3 4.40 3.20
## 4 3.98 3.13
## 5 4.94 3.53
## 6 3.18 2.23
## 7 4.15 2.60
## 8 3.66 2.71
## 9 4.23 3.06
## 10 4.80 3.59menghitung rata-rata 10 folds
mean_metric_polyA2 <- colMeans(metric_polyA2)
mean_metric_polyA2## rmse mae
## 4.323459 3.160102- Regresi polinomial berderajat 3
library(rsample)
set.seed(12345)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyA3 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(displacement,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_polyA3## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.14 3.14
## 2 5.40 4.14
## 3 3.39 2.53
## 4 6.05 3.95
## 5 4.05 3.07
## 6 3.35 2.61
## 7 5.10 3.33
## 8 4.39 3.67
## 9 3.04 2.43
## 10 3.66 2.83menghitung rata-rata 10 folds
mean_metric_polyA3 <- colMeans(metric_polyA3)
mean_metric_polyA3## rmse mae
## 4.257063 3.170277- Regresi polinomial berderajat 4
library(rsample)
set.seed(123456)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyA4 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(displacement,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_polyA4## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.53 3.47
## 2 5.83 4.35
## 3 3.76 2.54
## 4 4.37 2.79
## 5 4.23 3.16
## 6 3.86 2.89
## 7 5.18 3.89
## 8 3.16 2.35
## 9 4.01 3.12
## 10 4.26 3.19menghitung rata-rata 10 folds
mean_metric_polyA4 <- colMeans(metric_polyA4)
mean_metric_polyA4## rmse mae
## 4.318202 3.174789- Regresi polinomial berderajat 5
library(rsample)
set.seed(1234567)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyA5 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(displacement,5,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_polyA5## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 5.35 3.88
## 2 4.21 3.16
## 3 4.34 3.54
## 4 3.75 2.72
## 5 5.20 3.50
## 6 5.50 3.36
## 7 3.99 2.99
## 8 4.26 3.44
## 9 3.27 2.48
## 10 3.56 2.84menghitung rata-rata 10 folds
mean_metric_polyA5 <- colMeans(metric_polyA5)
mean_metric_polyA5## rmse mae
## 4.342347 3.189694- komparasi model regresi polinomial terbaik
## polynom rmse mae
## poly 1 4.562513 3.511871
## poly 2 4.323459 3.160102
## poly 3 4.257063 3.170277
## poly 4 4.318202 3.174789
## poly 5 4.342347 3.189694- model berdasarkan rmse
polynomA %>% slice_min(rmse)## polynom rmse mae
## mean_metric_polyA3 poly 3 4.257063 3.170277- model berdasarkan mae
polynomA %>% slice_min(mae)## polynom rmse mae
## mean_metric_polyA2 poly 2 4.323459 3.160102diperoleh model regresi polinomial berderajat 3 dengan rmse terkecil dan polinomial berderajat 2 dengan mae terkecil
Regresi fungsi tangga (visual)
- breaks = 2
model_tanggaA2 = lm(mpg ~ cut(displacement,2),data=Auto)
summary(model_tanggaA2)##
## Call:
## lm(formula = mpg ~ cut(displacement, 2), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3656 -4.1406 -0.4496 3.4833 22.9663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.3656 0.3536 74.56 <2e-16 ***
## cut(displacement, 2)(262,455] -11.3320 0.6966 -16.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.032 on 390 degrees of freedom
## Multiple R-squared: 0.4042, Adjusted R-squared: 0.4027
## F-statistic: 264.6 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,2),
lty = 1, col = "orange",se = F)+
theme_bw()
* breaks = 3 model_tanggaA3 = lm(mpg ~ cut(displacement,3),data=Auto)
summary(model_tanggaA3)##
## Call:
## lm(formula = mpg ~ cut(displacement, 3), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.4423 -2.8351 -0.2679 2.3649 20.1649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.6423 0.3320 86.27 <2e-16 ***
## cut(displacement, 3)(197,326] -10.8073 0.5700 -18.96 <2e-16 ***
## cut(displacement, 3)(326,455] -14.3745 0.7397 -19.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.947 on 389 degrees of freedom
## Multiple R-squared: 0.6004, Adjusted R-squared: 0.5983
## F-statistic: 292.2 on 2 and 389 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,3),
lty = 1, col = "orange",se = F)+
theme_bw()
* breaks = 4model_tanggaA4 = lm(mpg ~ cut(displacement,4),data=Auto)
summary(model_tanggaA4)##
## Call:
## lm(formula = mpg ~ cut(displacement, 4), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6099 -2.8099 -0.6897 2.3093 22.4250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.8099 0.3367 85.56 <2e-16 ***
## cut(displacement, 4)(165,262] -9.1188 0.6504 -14.02 <2e-16 ***
## cut(displacement, 4)(262,358] -13.2349 0.6699 -19.76 <2e-16 ***
## cut(displacement, 4)(358,455] -15.1202 0.9727 -15.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.914 on 388 degrees of freedom
## Multiple R-squared: 0.6066, Adjusted R-squared: 0.6036
## F-statistic: 199.4 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,4),
lty = 1, col = "orange",se = F)+
theme_bw()
* breaks = 5 model_tanggaA5 = lm(mpg ~ cut(displacement,5),data=Auto)
summary(model_tanggaA5)##
## Call:
## lm(formula = mpg ~ cut(displacement, 5), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2766 -2.8016 -0.2766 2.4493 18.9917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.2766 0.3400 86.104 < 2e-16 ***
## cut(displacement, 5)(145,223] -5.7528 0.8026 -7.168 3.88e-12 ***
## cut(displacement, 5)(223,300] -10.2682 0.6968 -14.736 < 2e-16 ***
## cut(displacement, 5)(300,378] -14.2259 0.6478 -21.959 < 2e-16 ***
## cut(displacement, 5)(378,455] -15.5766 1.0017 -15.549 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.711 on 387 degrees of freedom
## Multiple R-squared: 0.6394, Adjusted R-squared: 0.6356
## F-statistic: 171.5 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "orange",se = F)+
theme_bw()
* breaks = 6 model_tanggaA6 = lm(mpg ~ cut(displacement,6),data=Auto)
summary(model_tanggaA2)##
## Call:
## lm(formula = mpg ~ cut(displacement, 2), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3656 -4.1406 -0.4496 3.4833 22.9663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.3656 0.3536 74.56 <2e-16 ***
## cut(displacement, 2)(262,455] -11.3320 0.6966 -16.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.032 on 390 degrees of freedom
## Multiple R-squared: 0.4042, Adjusted R-squared: 0.4027
## F-statistic: 264.6 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,6),
lty = 1, col = "orange",se = F)+
theme_bw()
Regresi fungsi tangga (empiris)
set.seed(12345678)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
breaks <- 2:9
best_tanggaP1 <- map_dfr(breaks, function(i){
metric_tanggaP1 <- map_dfr(cross_val$splits,
function(x){
training <- Auto[x$in_id,]
training$displacement <- cut(training$displacement,i)
mod <- lm(mpg ~ displacement,
data=training)
labs_x <- levels(mod$model[,2])
labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))
testing <- Auto[-x$in_id,]
displacement_new <- cut(testing$displacement,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
pred <- predict(mod,
newdata=list(displacement=displacement_new))
truth <- testing$mpg
data_eval <- na.omit(data.frame(truth,pred))
rmse <- mlr3measures::rmse(truth = data_eval$truth,
response = data_eval$pred
)
mae <- mlr3measures::mae(truth = data_eval$truth,
response = data_eval$pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_tanggaP1
# menghitung rata-rata untuk 10 folds
mean_metric_tangga1 <- colMeans(metric_tanggaP1)
mean_metric_tangga1
})
best_tanggaP1 <- cbind(breaks=breaks,best_tanggaP1)
# menampilkan hasil all breaks
best_tanggaP1## breaks rmse mae
## 1 2 5.909043 4.627321
## 2 3 4.907504 3.683336
## 3 4 4.824395 3.623478
## 4 5 4.652961 3.488932
## 5 6 4.592822 3.398382
## 6 7 4.561135 3.366569
## 7 8 4.386949 3.239000
## 8 9 4.215282 3.090526- Rata-rata untuk 10 folds
mean_metric_tangga1<-colMeans(best_tanggaP1)
mean_metric_tangga1<-mean_metric_tangga1[-1]
mean_metric_tangga1## rmse mae
## 4.756261 3.564693#berdasarkan rmse
best_tanggaP1 %>% slice_min(rmse)## breaks rmse mae
## 1 9 4.215282 3.090526#berdasarkan mae
best_tanggaP1 %>% slice_min(mae)## breaks rmse mae
## 1 9 4.215282 3.090526diperoleh regresi fungsi tangga terbaik dengan breaks = 9
ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,9),
lty = 1, col = "orange",se = F)+
theme_bw()
Regresi natural cubic spline (visual)
* nilai knots yang ditentukan komputer (df=3)attr(ns(Auto$displacement, df=3),"knots")## 33.33333% 66.66667%
## 119 232library(splines)
model_splinedis3 = lm(mpg ~ ns(displacement, knots = c(119,232)),data=Auto)
summary(model_splinedis3)##
## Call:
## lm(formula = mpg ~ ns(displacement, knots = c(119, 232)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.8094 -2.3865 -0.2542 2.0600 20.3637
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.0671 0.7522 45.29 <2e-16 ***
## ns(displacement, knots = c(119, 232))1 -15.1851 1.0227 -14.85 <2e-16 ***
## ns(displacement, knots = c(119, 232))2 -27.7714 1.9402 -14.31 <2e-16 ***
## ns(displacement, knots = c(119, 232))3 -16.6492 1.1048 -15.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.36 on 388 degrees of freedom
## Multiple R-squared: 0.6903, Adjusted R-squared: 0.6879
## F-statistic: 288.3 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(119,232)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=4)attr(ns(Auto$displacement, df=4),"knots")## 25% 50% 75%
## 105.00 151.00 275.75library(splines)
model_splinedis4 = lm(mpg ~ ns(displacement, knots = c(105,151,275.75)),data=Auto)
summary(model_splinedis4)##
## Call:
## lm(formula = mpg ~ ns(displacement, knots = c(105, 151, 275.75)),
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.7432 -2.3595 -0.2352 2.0952 20.4131
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 33.999 1.070 31.761
## ns(displacement, knots = c(105, 151, 275.75))1 -12.157 1.268 -9.588
## ns(displacement, knots = c(105, 151, 275.75))2 -17.094 1.346 -12.696
## ns(displacement, knots = c(105, 151, 275.75))3 -24.535 2.502 -9.808
## ns(displacement, knots = c(105, 151, 275.75))4 -18.385 1.341 -13.706
## Pr(>|t|)
## (Intercept) <2e-16 ***
## ns(displacement, knots = c(105, 151, 275.75))1 <2e-16 ***
## ns(displacement, knots = c(105, 151, 275.75))2 <2e-16 ***
## ns(displacement, knots = c(105, 151, 275.75))3 <2e-16 ***
## ns(displacement, knots = c(105, 151, 275.75))4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.368 on 387 degrees of freedom
## Multiple R-squared: 0.6901, Adjusted R-squared: 0.6869
## F-statistic: 215.4 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(105,151,275.75)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=5)attr(ns(Auto$displacement, df=5),"knots")## 20% 40% 60% 80%
## 98 122 225 305library(splines)
model_splinedis5 = lm(mpg ~ ns(displacement, knots = c(98,122,225,305)),data=Auto)
summary(model_splinedis5)##
## Call:
## lm(formula = mpg ~ ns(displacement, knots = c(98, 122, 225, 305)),
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.4363 -2.3774 -0.3442 2.2338 19.9381
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 31.374 1.334 23.512
## ns(displacement, knots = c(98, 122, 225, 305))1 -9.512 1.447 -6.575
## ns(displacement, knots = c(98, 122, 225, 305))2 -10.204 1.819 -5.610
## ns(displacement, knots = c(98, 122, 225, 305))3 -18.122 1.299 -13.951
## ns(displacement, knots = c(98, 122, 225, 305))4 -16.960 2.957 -5.736
## ns(displacement, knots = c(98, 122, 225, 305))5 -18.097 1.457 -12.420
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## ns(displacement, knots = c(98, 122, 225, 305))1 1.59e-10 ***
## ns(displacement, knots = c(98, 122, 225, 305))2 3.87e-08 ***
## ns(displacement, knots = c(98, 122, 225, 305))3 < 2e-16 ***
## ns(displacement, knots = c(98, 122, 225, 305))4 1.96e-08 ***
## ns(displacement, knots = c(98, 122, 225, 305))5 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.326 on 386 degrees of freedom
## Multiple R-squared: 0.6967, Adjusted R-squared: 0.6928
## F-statistic: 177.3 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(98,122,225,305)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=6)attr(ns(Auto$displacement, df=6),"knots")## 16.66667% 33.33333% 50% 66.66667% 83.33333%
## 97 119 151 232 318library(splines)
model_splinedis6 = lm(mpg ~ ns(displacement, knots = c(97,119,151,232,318)),data=Auto)
summary(model_splinedis6)##
## Call:
## lm(formula = mpg ~ ns(displacement, knots = c(97, 119, 151, 232,
## 318)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.9381 -2.3829 -0.3102 2.1258 20.6988
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 29.254 1.409
## ns(displacement, knots = c(97, 119, 151, 232, 318))1 -5.603 1.473
## ns(displacement, knots = c(97, 119, 151, 232, 318))2 -2.879 2.154
## ns(displacement, knots = c(97, 119, 151, 232, 318))3 -11.521 1.735
## ns(displacement, knots = c(97, 119, 151, 232, 318))4 -16.379 1.372
## ns(displacement, knots = c(97, 119, 151, 232, 318))5 -10.755 3.232
## ns(displacement, knots = c(97, 119, 151, 232, 318))6 -19.107 1.483
## t value Pr(>|t|)
## (Intercept) 20.760 < 2e-16 ***
## ns(displacement, knots = c(97, 119, 151, 232, 318))1 -3.803 0.000166 ***
## ns(displacement, knots = c(97, 119, 151, 232, 318))2 -1.336 0.182237
## ns(displacement, knots = c(97, 119, 151, 232, 318))3 -6.639 1.08e-10 ***
## ns(displacement, knots = c(97, 119, 151, 232, 318))4 -11.935 < 2e-16 ***
## ns(displacement, knots = c(97, 119, 151, 232, 318))5 -3.327 0.000961 ***
## ns(displacement, knots = c(97, 119, 151, 232, 318))6 -12.885 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.24 on 385 degrees of freedom
## Multiple R-squared: 0.7094, Adjusted R-squared: 0.7049
## F-statistic: 156.6 on 6 and 385 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(97,119,151,232,318)),
lty = 1,se=F)
Regresi natural cubic spline (empiris)
library(splines)
set.seed(123456789)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 2:10
best_splineA <- map_dfr(df, function(i){
metric_splineA <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(displacement,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_splineA
# menghitung rata-rata untuk 10 folds
mean_metric_splineA <- colMeans(metric_splineA)
mean_metric_splineA
}
)best_splineA <- cbind(df=df,best_splineA)
best_splineA## df rmse mae
## 1 2 4.323392 3.171205
## 2 3 4.348383 3.191347
## 3 4 4.363702 3.203273
## 4 5 4.320739 3.212415
## 5 6 4.239470 3.142298
## 6 7 4.206470 3.134108
## 7 8 4.206469 3.154002
## 8 9 4.187126 3.153311
## 9 10 4.188307 3.154775mean_metric_splineA<-colMeans(best_splineA[-1])
mean_metric_splineA## rmse mae
## 4.264895 3.168526best_splineA %>% slice_min(rmse)## df rmse mae
## 1 9 4.187126 3.153311best_splineA %>% slice_min(mae)## df rmse mae
## 1 7 4.20647 3.134108model regresi natural cubic spline dengan df = 9 memiliki rmse terkecil dan model regresi natural cubic spline dengan df = 7 memiliki mae terkecil.
attr(ns(Auto$displacement, df=7),"knots")## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429%
## 96.85714 108.00000 134.57143 198.00000 250.00000 321.14286attr(ns(Auto$displacement, df=9),"knots")## 11.11111% 22.22222% 33.33333% 44.44444% 55.55556% 66.66667% 77.77778% 88.88889%
## 91.0000 98.0000 119.0000 140.0000 173.0000 232.0000 302.2222 350.0000
Komparasi model
- Natural cubic spline dengan knots
attr(ns(Auto$displacement, df=9),"knots")## 11.11111% 22.22222% 33.33333% 44.44444% 55.55556% 66.66667% 77.77778% 88.88889%
## 91.0000 98.0000 119.0000 140.0000 173.0000 232.0000 302.2222 350.0000set.seed(0987)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_splineA9 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(displacement, knots = c(91,98,119,140,173,232,302.2222,350)),
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_splineA9## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.06 2.41
## 2 4.53 3.55
## 3 4.85 3.77
## 4 3.26 2.64
## 5 4.93 3.89
## 6 4.07 2.92
## 7 4.45 3.34
## 8 4.86 2.97
## 9 4.33 3.18
## 10 3.34 2.70menghitung rata-rata 10 folds
mean_metric_splineA9 <- colMeans(metric_splineA9)
mean_metric_splineA9## rmse mae
## 4.166547 3.138199## model_displacement rmse mae
## Linear 4.575194 3.509608
## Poly 1 4.562513 3.511871
## Poly 2 4.323459 3.160102
## Poly 3 4.257063 3.170277
## Poly 4 4.318202 3.174789
## Poly 5 4.342347 3.189694
## Tangga 4.756261 3.564693
## Spline df =7 4.264895 3.168526
## Spline knots = 8 4.166547 3.138199BestAfinal %>% slice_min(rmse)## model_displacement rmse mae
## mean_metric_splineA9 Spline knots = 8 4.166547 3.138199BestAfinal %>% slice_min(mae)## model_displacement rmse mae
## mean_metric_splineA9 Spline knots = 8 4.166547 3.138199diperoleh model regresi untuk respons mpg dan prediktor displacementdengan nilai rmse dan mae terkecil yaitu model regresi natural cubic spline dengan knots = 8.
Soal 2 : mpg - weight
## Regresi linear
model_linearB=lm(mpg~weight,data=Auto)
summary(model_linearB)##
## Call:
## lm(formula = mpg ~ weight, data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.9736 -2.7556 -0.3358 2.1379 16.5194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.216524 0.798673 57.87 <2e-16 ***
## weight -0.007647 0.000258 -29.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.333 on 390 degrees of freedom
## Multiple R-squared: 0.6926, Adjusted R-squared: 0.6918
## F-statistic: 878.8 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.50, color="grey") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "pink",se = F)+
theme_bw() 
set.seed(09876)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_linearB <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ weight,
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_linearB## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.21 3.46
## 2 4.56 3.56
## 3 4.66 3.41
## 4 4.08 2.93
## 5 4.48 3.38
## 6 4.72 3.48
## 7 3.64 3.02
## 8 4.17 3.48
## 9 5.35 3.76
## 10 3.09 2.39mean_metric_linearB <- colMeans(metric_linearB)
mean_metric_linearB## rmse mae
## 4.295467 3.287539Regresi polinomial (visual)
- Regresi polinomial derajat 1
model_polinomialB = lm(mpg~poly(weight,1,raw = T),data=Auto)
summary(model_polinomialB)##
## Call:
## lm(formula = mpg ~ poly(weight, 1, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.9736 -2.7556 -0.3358 2.1379 16.5194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.216524 0.798673 57.87 <2e-16 ***
## poly(weight, 1, raw = T) -0.007647 0.000258 -29.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.333 on 390 degrees of freedom
## Multiple R-squared: 0.6926, Adjusted R-squared: 0.6918
## F-statistic: 878.8 on 1 and 390 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,1,raw=T),
lty = 1, col = "pink",se = F)+
theme_bw()
- Regresi polinomial derajat 2
model_polinomialB2 = lm(mpg~poly(weight,2,raw = T),data=Auto)
summary(model_polinomialB2)##
## Call:
## lm(formula = mpg ~ poly(weight, 2, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6246 -2.7134 -0.3485 1.8267 16.0866
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.226e+01 2.993e+00 20.800 < 2e-16 ***
## poly(weight, 2, raw = T)1 -1.850e-02 1.972e-03 -9.379 < 2e-16 ***
## poly(weight, 2, raw = T)2 1.697e-06 3.059e-07 5.545 5.43e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.176 on 389 degrees of freedom
## Multiple R-squared: 0.7151, Adjusted R-squared: 0.7137
## F-statistic: 488.3 on 2 and 389 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "pink",se = F)+
theme_bw()
- Regresi polinomial derajat 3
model_polinomialB3 = lm(mpg ~ poly(weight,3,raw = T),data=Auto)
summary(model_polinomialB3)##
## Call:
## lm(formula = mpg ~ poly(weight, 3, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6259 -2.7080 -0.3552 1.8385 16.0816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.170e+01 1.104e+01 5.587 4.36e-08 ***
## poly(weight, 3, raw = T)1 -1.793e-02 1.091e-02 -1.643 0.101
## poly(weight, 3, raw = T)2 1.515e-06 3.450e-06 0.439 0.661
## poly(weight, 3, raw = T)3 1.846e-11 3.503e-10 0.053 0.958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.182 on 388 degrees of freedom
## Multiple R-squared: 0.7151, Adjusted R-squared: 0.7129
## F-statistic: 324.7 on 3 and 388 DF, p-value: < 2.2e-16 ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,3,raw=T),
lty = 1, col = "pink",se = F)+
theme_bw()
- Regresi polinomial derajat 4
model_polinomialB4 = lm(mpg~poly(weight,4,raw = T),data=Auto)
summary(model_polinomialB4)##
## Call:
## lm(formula = mpg ~ poly(weight, 4, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7233 -2.7179 -0.4254 1.9112 16.1935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.685e+01 3.864e+01 0.954 0.341
## poly(weight, 4, raw = T)1 1.548e-02 5.097e-02 0.304 0.762
## poly(weight, 4, raw = T)2 -1.470e-05 2.441e-05 -0.602 0.547
## poly(weight, 4, raw = T)3 3.389e-09 5.035e-09 0.673 0.501
## poly(weight, 4, raw = T)4 -2.539e-13 3.784e-13 -0.671 0.503
##
## Residual standard error: 4.185 on 387 degrees of freedom
## Multiple R-squared: 0.7155, Adjusted R-squared: 0.7125
## F-statistic: 243.3 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "pink",se = F)+
theme_bw()
- Regresi polinomial derajat 5
model_polinomialB5 = lm(mpg~poly(weight,5,raw = T),data=Auto)
summary(model_polinomialB5)##
## Call:
## lm(formula = mpg ~ poly(weight, 5, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.8872 -2.7327 -0.3687 1.8744 16.3456
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.083e+01 1.343e+02 -0.602 0.548
## poly(weight, 5, raw = T)1 2.149e-01 2.238e-01 0.960 0.338
## poly(weight, 5, raw = T)2 -1.456e-04 1.451e-04 -1.003 0.316
## poly(weight, 5, raw = T)3 4.505e-08 4.580e-08 0.984 0.326
## poly(weight, 5, raw = T)4 -6.685e-12 7.038e-12 -0.950 0.343
## poly(weight, 5, raw = T)5 3.861e-16 4.220e-16 0.915 0.361
##
## Residual standard error: 4.186 on 386 degrees of freedom
## Multiple R-squared: 0.7161, Adjusted R-squared: 0.7124
## F-statistic: 194.7 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,5,raw=T),
lty = 1, col = "pink",se = F)+
theme_bw()
Regresi polinomial (empiris)
- Regresi polinomial berderajat 1
library(rsample)
set.seed(098765)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyB1 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(weight,1,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_polyB1## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.85 3.19
## 2 4.83 3.46
## 3 4.82 3.54
## 4 4.11 2.75
## 5 4.29 3.49
## 6 4.06 3.17
## 7 4.41 3.55
## 8 4.44 3.53
## 9 4.47 3.35
## 10 3.89 2.83menghitung rata-rata 10 folds
mean_metric_polyB1 <- colMeans(metric_polyB1)
mean_metric_polyB1## rmse mae
## 4.315989 3.285392- Regresi polinomial berderajat 2
library(rsample)
set.seed(087654)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyB2 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(weight,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_polyB2## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.85 2.99
## 2 3.95 2.91
## 3 3.80 2.97
## 4 4.57 3.17
## 5 3.70 2.89
## 6 4.22 3.12
## 7 3.66 2.80
## 8 4.23 3.16
## 9 5.08 3.62
## 10 4.59 3.10menghitung rata-rata 10 folds
mean_metric_polyB2 <- colMeans(metric_polyB2)
mean_metric_polyB2## rmse mae
## 4.164003 3.071097- Regresi polinomial berderajat 3
library(rsample)
set.seed(09876543)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyB3 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(weight,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_polyB3## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.42 3.37
## 2 3.38 2.56
## 3 5.07 3.39
## 4 4.05 2.95
## 5 3.74 2.83
## 6 4.83 3.36
## 7 4.83 3.55
## 8 3.35 2.66
## 9 3.56 2.85
## 10 4.16 3.23menghitung rata-rata 10 folds
mean_metric_polyB3 <- colMeans(metric_polyB3)
mean_metric_polyB3## rmse mae
## 4.140591 3.075079- Regresi polinomial berderajat 4
library(rsample)
set.seed(098765432)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyB4 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(weight,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_polyB4## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.41 3.23
## 2 4.87 3.92
## 3 5.07 3.48
## 4 3.83 2.98
## 5 2.76 2.10
## 6 3.02 2.49
## 7 5.23 4.08
## 8 3.63 2.65
## 9 3.11 2.30
## 10 5.14 3.57menghitung rata-rata 10 folds
mean_metric_polyB4 <- colMeans(metric_polyB4)
mean_metric_polyB4## rmse mae
## 4.106271 3.079978- Regresi polinomial berderajat 5
library(rsample)
set.seed(0987654321)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyB5 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(weight,5,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_polyB5## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.23 3.02
## 2 3.82 2.87
## 3 4.66 3.56
## 4 3.20 2.38
## 5 4.53 3.31
## 6 3.78 3.02
## 7 4.91 3.44
## 8 4.20 3.14
## 9 4.37 3.21
## 10 4.13 3.01menghitung rata-rata 10 folds
mean_metric_polyB5 <- colMeans(metric_polyB5)
mean_metric_polyB5## rmse mae
## 4.182899 3.094782- komparasi model regresi polinomial terbaik
## polynom rmse mae
## poly 1 4.315989 3.285392
## poly 2 4.164003 3.071097
## poly 3 4.140591 3.075079
## poly 4 4.106271 3.079978
## poly 5 4.182899 3.094782- model berdasarkan rmse
polynom_Wei %>% slice_min(rmse)## polynom rmse mae
## mean_metric_polyB4 poly 4 4.106271 3.079978- model berdasarkan mae
polynom_Wei %>% slice_min(mae)## polynom rmse mae
## mean_metric_polyB2 poly 2 4.164003 3.071097diperoleh model regresi polinomial berderajat 4 dengan rmse terkecil dan polinomial berderajat 2 dengan mae terkecil
ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, aes(col = "Derajat 2"),se = F)+
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, aes(col = "Derajat 4"),se = F)+labs(color="regresi polynomial")+
scale_color_manual(values = c("Derajat 2"="red","Derajat 4"="blue"))+theme_bw()
Regresi fungsi tangga (visual)
- breaks = 2
model_tanggaB2 = lm(mpg ~ cut(weight,2),data=Auto)
summary(model_tanggaB2)##
## Call:
## lm(formula = mpg ~ cut(weight, 2), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1913 -3.3163 -0.6307 3.3087 19.4087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.1913 0.3456 78.68 <2e-16 ***
## cut(weight, 2)(3.38e+03,5.14e+03] -11.5606 0.6072 -19.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.626 on 390 degrees of freedom
## Multiple R-squared: 0.4817, Adjusted R-squared: 0.4804
## F-statistic: 362.5 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,2),
lty = 1, col = "purple",se = F)+
theme_bw()
* breaks = 3 model_tanggaB3 = lm(mpg ~ cut(weight,3),data=Auto)
summary(model_tanggaB3)##
## Call:
## lm(formula = mpg ~ cut(weight, 3), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3604 -3.3604 -0.4422 2.3578 18.3578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.3604 0.3509 83.67 <2e-16 ***
## cut(weight, 3)(2.79e+03,3.96e+03] -9.7182 0.5461 -17.80 <2e-16 ***
## cut(weight, 3)(3.96e+03,5.14e+03] -15.4850 0.6977 -22.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.862 on 389 degrees of freedom
## Multiple R-squared: 0.6139, Adjusted R-squared: 0.6119
## F-statistic: 309.3 on 2 and 389 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,3),
lty = 1, col = "purple",se = F)+
theme_bw()
* breaks = 4model_tanggaB4 = lm(mpg ~ cut(weight,4),data=Auto)
summary(model_tanggaB4)##
## Call:
## lm(formula = mpg ~ cut(weight, 4), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.541 -3.065 -0.339 2.477 16.059
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.5413 0.3923 77.84 <2e-16 ***
## cut(weight, 4)(2.49e+03,3.38e+03] -7.2765 0.5782 -12.58 <2e-16 ***
## cut(weight, 4)(3.38e+03,4.26e+03] -13.8180 0.6402 -21.58 <2e-16 ***
## cut(weight, 4)(4.26e+03,5.14e+03] -17.2022 0.8312 -20.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.692 on 388 degrees of freedom
## Multiple R-squared: 0.6414, Adjusted R-squared: 0.6386
## F-statistic: 231.3 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,4),
lty = 1, col = "purple",se = F)+
theme_bw()
* breaks = 5 model_tanggaB5 = lm(mpg ~ cut(displacement,5),data=Auto)
summary(model_tanggaA5)##
## Call:
## lm(formula = mpg ~ cut(displacement, 5), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2766 -2.8016 -0.2766 2.4493 18.9917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.2766 0.3400 86.104 < 2e-16 ***
## cut(displacement, 5)(145,223] -5.7528 0.8026 -7.168 3.88e-12 ***
## cut(displacement, 5)(223,300] -10.2682 0.6968 -14.736 < 2e-16 ***
## cut(displacement, 5)(300,378] -14.2259 0.6478 -21.959 < 2e-16 ***
## cut(displacement, 5)(378,455] -15.5766 1.0017 -15.549 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.711 on 387 degrees of freedom
## Multiple R-squared: 0.6394, Adjusted R-squared: 0.6356
## F-statistic: 171.5 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=displacement, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "orange",se = F)+
theme_bw()
* breaks = 6 model_tanggaB6 = lm(mpg ~ cut(weight,6),data=Auto)
summary(model_tanggaB6)##
## Call:
## lm(formula = mpg ~ cut(weight, 6), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.603 -2.486 -0.423 1.635 16.901
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.6033 0.4433 73.55 <2e-16 ***
## cut(weight, 6)(2.2e+03,2.79e+03] -6.1043 0.6082 -10.04 <2e-16 ***
## cut(weight, 6)(2.79e+03,3.38e+03] -11.1170 0.6624 -16.78 <2e-16 ***
## cut(weight, 6)(3.38e+03,3.96e+03] -15.1324 0.6940 -21.80 <2e-16 ***
## cut(weight, 6)(3.96e+03,4.55e+03] -18.2380 0.7466 -24.43 <2e-16 ***
## cut(weight, 6)(4.55e+03,5.14e+03] -20.2283 1.1409 -17.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.205 on 386 degrees of freedom
## Multiple R-squared: 0.7134, Adjusted R-squared: 0.7097
## F-statistic: 192.2 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,6),
lty = 1, col = "purple",se = F)+
theme_bw()
Regresi fungsi tangga (empiris)
set.seed(112233)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
breaks <- 2:10
best_tanggaP2 <- map_dfr(breaks, function(i){
metric_tanggaP2 <- map_dfr(cross_val$splits,
function(x){
training <- Auto[x$in_id,]
training$weight <- cut(training$weight,i)
mod <- lm(mpg ~ weight,
data=training)
labs_x <- levels(mod$model[,2])
labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))
testing <- Auto[-x$in_id,]
weight_new <- cut(testing$weight,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
pred <- predict(mod,
newdata=list(weight=weight_new))
truth <- testing$mpg
data_eval <- na.omit(data.frame(truth,pred))
rmse <- mlr3measures::rmse(truth = data_eval$truth,
response = data_eval$pred
)
mae <- mlr3measures::mae(truth = data_eval$truth,
response = data_eval$pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_tanggaP2
# menghitung rata-rata untuk 10 folds
mean_metric_tangga2 <- colMeans(metric_tanggaP2)
mean_metric_tangga2
})
best_tanggaP2 <- cbind(breaks=breaks,best_tanggaP2)
# menampilkan hasil all breaks
best_tanggaP2## breaks rmse mae
## 1 2 5.606718 4.360467
## 2 3 4.845242 3.654466
## 3 4 4.669452 3.602658
## 4 5 4.454064 3.259862
## 5 6 4.213537 3.109653
## 6 7 4.313684 3.230943
## 7 8 4.360573 3.237396
## 8 9 4.285401 3.142870
## 9 10 4.354320 3.188635- Rata-rata untuk 10 folds
mean_metric_tangga2<-colMeans(best_tanggaP2)
mean_metric_tangga2<-mean_metric_tangga2[-1]
mean_metric_tangga2## rmse mae
## 4.566999 3.420772#berdasarkan rmse
best_tanggaP2 %>% slice_min(rmse)## breaks rmse mae
## 1 6 4.213537 3.109653#berdasarkan mae
best_tanggaP2 %>% slice_min(mae)## breaks rmse mae
## 1 6 4.213537 3.109653diperoleh regresi fungsi tangga terbaik dengan breaks = 6
ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,6),
lty = 1, col = "purple",se = F)+
theme_bw()
Regresi natural cubic spline (visual)
* nilai knots yang ditentukan komputer (df=3)attr(ns(Auto$weight, df=3),"knots")## 33.33333% 66.66667%
## 2397.000 3333.667library(splines)
model_splinewei3 = lm(mpg ~ ns(weight, knots = c(2397,3333.667)),data=Auto)
summary(model_splinewei3)##
## Call:
## lm(formula = mpg ~ ns(weight, knots = c(2397, 3333.667)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6780 -2.6926 -0.4112 1.8169 16.1671
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.5397 0.9896 36.92 <2e-16 ***
## ns(weight, knots = c(2397, 3333.667))1 -17.2157 0.9686 -17.77 <2e-16 ***
## ns(weight, knots = c(2397, 3333.667))2 -31.5415 2.3954 -13.17 <2e-16 ***
## ns(weight, knots = c(2397, 3333.667))3 -20.1636 1.2138 -16.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.18 on 388 degrees of freedom
## Multiple R-squared: 0.7154, Adjusted R-squared: 0.7132
## F-statistic: 325.1 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(2397,3333.667)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=4)attr(ns(Auto$weight, df=4),"knots")## 25% 50% 75%
## 2225.25 2803.50 3614.75library(splines)
model_splinewei4 = lm(mpg ~ ns(weight, knots = c(2225.25,2803.50,3614.75)),data=Auto)
summary(model_splinewei4)##
## Call:
## lm(formula = mpg ~ ns(weight, knots = c(2225.25, 2803.5, 3614.75)),
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7082 -2.7027 -0.3953 1.8192 16.1999
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 36.377 1.354 26.863
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))1 -13.793 1.314 -10.494
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))2 -18.746 1.348 -13.902
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))3 -29.156 3.126 -9.326
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))4 -21.514 1.488 -14.456
## Pr(>|t|)
## (Intercept) <2e-16 ***
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))1 <2e-16 ***
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))2 <2e-16 ***
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))3 <2e-16 ***
## ns(weight, knots = c(2225.25, 2803.5, 3614.75))4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.185 on 387 degrees of freedom
## Multiple R-squared: 0.7154, Adjusted R-squared: 0.7125
## F-statistic: 243.2 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(2225.25,2803.50,3614.75)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=5)attr(ns(Auto$weight, df=5),"knots")## 20% 40% 60% 80%
## 2155.0 2583.2 3113.4 3820.8library(splines)
model_splinewei5 = lm(mpg ~ ns(weight, knots = c(2155,2583.2,3113.4,3820.8)),data=Auto)
summary(model_splinewei5)##
## Call:
## lm(formula = mpg ~ ns(weight, knots = c(2155, 2583.2, 3113.4,
## 3820.8)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.8856 -2.7740 -0.4045 1.9150 16.1559
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 35.686 1.591
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))1 -10.885 1.553
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))2 -15.031 1.991
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))3 -19.698 1.433
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))4 -26.677 3.648
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))5 -21.759 1.674
## t value Pr(>|t|)
## (Intercept) 22.423 < 2e-16 ***
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))1 -7.010 1.07e-11 ***
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))2 -7.549 3.20e-13 ***
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))3 -13.745 < 2e-16 ***
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))4 -7.312 1.53e-12 ***
## ns(weight, knots = c(2155, 2583.2, 3113.4, 3820.8))5 -12.997 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.188 on 386 degrees of freedom
## Multiple R-squared: 0.7158, Adjusted R-squared: 0.7121
## F-statistic: 194.4 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(2155,2583.2,3113.4,3820.8)),
lty = 1,se=F)
* nilai knots yang ditentukan komputer (df=6)attr(ns(Auto$weight, df=6),"knots")## 16.66667% 33.33333% 50% 66.66667% 83.33333%
## 2125.000 2397.000 2803.500 3333.667 3960.833library(splines)
model_splinewei6 = lm(mpg ~ ns(weight, knots = c(2125,2397,2803.5,3333.667,3960.833)),data=Auto)
summary(model_splinewei6)##
## Call:
## lm(formula = mpg ~ ns(weight, knots = c(2125, 2397, 2803.5, 3333.667,
## 3960.833)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.318 -2.661 -0.450 1.964 15.637
##
## Coefficients:
## Estimate
## (Intercept) 33.697
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))1 -8.602
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))2 -8.915
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))3 -15.623
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))4 -18.733
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))5 -21.079
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))6 -22.688
## Std. Error
## (Intercept) 1.811
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))1 1.754
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))2 2.295
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))3 2.001
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))4 1.602
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))5 4.120
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))6 1.843
## t value Pr(>|t|)
## (Intercept) 18.602 < 2e-16
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))1 -4.905 1.38e-06
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))2 -3.884 0.000121
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))3 -7.806 5.62e-14
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))4 -11.697 < 2e-16
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))5 -5.116 4.92e-07
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))6 -12.312 < 2e-16
##
## (Intercept) ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))1 ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))2 ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))3 ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))4 ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))5 ***
## ns(weight, knots = c(2125, 2397, 2803.5, 3333.667, 3960.833))6 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.167 on 385 degrees of freedom
## Multiple R-squared: 0.7193, Adjusted R-squared: 0.715
## F-statistic: 164.5 on 6 and 385 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=weight, y=mpg)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(2125,2397,2803.5,3333.667,3960.833)),
lty = 1,se=F)
Regresi natural cubic spline (empiris)
library(splines)
set.seed(11223344)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 2:10
best_splineB <- map_dfr(df, function(i){
metric_splineB <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(weight,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_splineB
# menghitung rata-rata untuk 10 folds
mean_metric_splineB <- colMeans(metric_splineB)
mean_metric_splineB
}
)best_splineB <- cbind(df=df,best_splineB)
best_splineB## df rmse mae
## 1 2 4.110309 3.060302
## 2 3 4.116047 3.063867
## 3 4 4.124727 3.069173
## 4 5 4.128629 3.080111
## 5 6 4.116874 3.085895
## 6 7 4.100994 3.045973
## 7 8 4.114049 3.064631
## 8 9 4.114405 3.054340
## 9 10 4.112056 3.042977mean_metric_splineB<-colMeans(best_splineB[-1])
mean_metric_splineB## rmse mae
## 4.115343 3.063030best_splineB %>% slice_min(rmse)## df rmse mae
## 1 7 4.100994 3.045973best_splineB %>% slice_min(mae)## df rmse mae
## 1 10 4.112056 3.042977model regresi natural cubic spline dengan df = 7 memiliki rmse terkecil dan model regresi natural cubic spline dengan df = 10 memiliki mae terkecil.
attr(ns(Auto$weight, df=7),"knots")## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429%
## 2074.857 2285.429 2635.000 2986.571 3446.143 4096.286attr(ns(Auto$weight, df=10),"knots")## 10% 20% 30% 40% 50% 60% 70% 80% 90%
## 1990.0 2155.0 2303.0 2583.2 2803.5 3113.4 3428.5 3820.8 4277.6
Komparasi model
- Natural cubic spline dengan knots berdasarkan rmse terkecil
attr(ns(Auto$weight, df=7),"knots")## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429%
## 2074.857 2285.429 2635.000 2986.571 3446.143 4096.286set.seed(1122334455)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_splineB7 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(weight, knots = c(2074.857,2285.429,2635,2986.571,3446.143, 4096.286)),
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_splineB7## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.74 2.92
## 2 3.38 2.62
## 3 4.40 3.03
## 4 4.39 3.01
## 5 4.37 3.47
## 6 3.58 2.76
## 7 5.34 3.90
## 8 4.76 3.23
## 9 4.05 3.04
## 10 3.20 2.62menghitung rata-rata 10 folds
mean_metric_splineB7 <- colMeans(metric_splineB7)
mean_metric_splineB7## rmse mae
## 4.119871 3.060278## model_weight rmse mae
## Linear 4.295467 3.287539
## Poly 1 4.315989 3.285392
## Poly 2 4.164003 3.071097
## Poly 3 4.140591 3.075079
## Poly 4 4.106271 3.079978
## Poly 5 4.182899 3.094782
## Tangga 4.566999 3.420772
## Spline df=7 4.115343 3.063030
## Spline knots = 6 4.119871 3.060278BestBfinal %>% slice_min(rmse)## model_weight rmse mae
## mean_metric_polyB4 Poly 4 4.106271 3.079978BestBfinal %>% slice_min(mae)## model_weight rmse mae
## mean_metric_splineB7 Spline knots = 6 4.119871 3.060278diperoleh model regresi untuk respons mpg dan prediktor weightdengan nilai rmse terkecil adalah regresi polynomial berderajat 4 dan mae terkecil yaitu model regresi natural cubic spline dengan knots = 6.
Soal 2 : mpg - horsepower
Regresi linear
model_linearC=lm(mpg~horsepower,data=Auto)
summary(model_linearC)##
## 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-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.50, color="grey") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "red",se = F)+
theme_bw() 
Regresi linear (Empiris)
set.seed(980510)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_linearC <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ horsepower,
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_linearC## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 5.43 4.07
## 2 4.89 3.43
## 3 5.23 4.37
## 4 5.24 4.05
## 5 5.69 4.46
## 6 4.22 3.47
## 7 4.78 3.65
## 8 4.30 3.45
## 9 4.62 3.83
## 10 4.45 3.57Rata-rata untuk 10 folds
mean_metric_linearC <- colMeans(metric_linearC)
mean_metric_linearC## rmse mae
## 4.883887 3.835224Regresi Polinomial (Visual)
- Regresi polinomial derajat 1
model_polinomialC1 = lm(mpg~poly(horsepower,1,raw = T),data=Auto)
summary(model_polinomialC1)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 1, raw = T), 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 ***
## poly(horsepower, 1, raw = T) -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 ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,1,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
- Regresi polinomial derajat 2
model_polinomialC2 = lm(mpg~poly(horsepower,2,raw = T),data=Auto)
summary(model_polinomialC2)##
## 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 ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
- Regresi polinomial derajat 3
model_polinomialC3 = lm(mpg ~ poly(horsepower,3,raw = T),data=Auto)
summary(model_polinomialC3)##
## 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 ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,3,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
- Regresi polinomial derajat 4
model_polinomialC4 = lm(mpg~poly(horsepower,4,raw = T),data=Auto)
summary(model_polinomialC4)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 4, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.8820 -2.5802 -0.1682 2.2100 16.1434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.757e+01 1.196e+01 3.977 8.32e-05 ***
## poly(horsepower, 4, raw = T)1 -7.667e-02 4.313e-01 -0.178 0.859
## poly(horsepower, 4, raw = T)2 -4.345e-03 5.497e-03 -0.790 0.430
## poly(horsepower, 4, raw = T)3 3.245e-05 2.926e-05 1.109 0.268
## poly(horsepower, 4, raw = T)4 -6.530e-08 5.504e-08 -1.186 0.236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.373 on 387 degrees of freedom
## Multiple R-squared: 0.6893, Adjusted R-squared: 0.6861
## F-statistic: 214.7 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
- Regresi polinomial derajat 5
model_polinomialC5 = lm(mpg~poly(horsepower,5,raw = T),data=Auto)
summary(model_polinomialC5)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 5, raw = T), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.4326 -2.5285 -0.2925 2.1750 15.9730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.223e+01 2.857e+01 -1.128 0.26003
## poly(horsepower, 5, raw = T)1 3.700e+00 1.303e+00 2.840 0.00475 **
## poly(horsepower, 5, raw = T)2 -7.142e-02 2.253e-02 -3.170 0.00164 **
## poly(horsepower, 5, raw = T)3 5.931e-04 1.850e-04 3.206 0.00146 **
## poly(horsepower, 5, raw = T)4 -2.281e-06 7.243e-07 -3.150 0.00176 **
## poly(horsepower, 5, raw = T)5 3.330e-09 1.085e-09 3.068 0.00231 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.326 on 386 degrees of freedom
## Multiple R-squared: 0.6967, Adjusted R-squared: 0.6928
## F-statistic: 177.4 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,5,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
Regresi polinomial (Empiris)
- Regresi polinomial derajat 1
library(rsample)
set.seed(15)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyC1 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(horsepower,1,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_polyC1## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 5.56 4.21
## 2 4.45 3.49
## 3 4.17 3.33
## 4 5.61 4.12
## 5 4.66 3.76
## 6 5.54 4.44
## 7 4.65 3.49
## 8 4.49 3.58
## 9 5.19 4.27
## 10 4.45 3.67menghitung rata-rata 10 folds
mean_metric_polyC1 <- colMeans(metric_polyC1)
mean_metric_polyC1## rmse mae
## 4.875532 3.835955- Regresi polinomial derajat 2
library(rsample)
set.seed(150)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyC2 <- 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_polyC2## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.78 3.47
## 2 5.29 3.88
## 3 3.38 2.62
## 4 4.59 3.22
## 5 4.91 3.59
## 6 4.87 3.89
## 7 3.27 2.71
## 8 5.01 3.84
## 9 3.26 2.45
## 10 3.98 3.14menghitung rata-rata 10 folds
mean_metric_polyC2 <- colMeans(metric_polyC2)
mean_metric_polyC2## rmse mae
## 4.334441 3.280627- Regresi polinomial derajat 3
set.seed(15012)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyC3 <- 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_polyC3## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 6.33 4.75
## 2 4.95 3.61
## 3 3.71 3.01
## 4 3.32 2.44
## 5 4.18 2.95
## 6 4.85 3.73
## 7 3.81 3.06
## 8 3.50 2.68
## 9 3.78 2.96
## 10 4.43 3.41menghitung rata-rata untuk 10 folds
mean_metric_polyC3 <- colMeans(metric_polyC3)
mean_metric_polyC3## rmse mae
## 4.286842 3.259572- Regresi polinomial derajat 4
set.seed(150121)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyC4 <- 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_polyC4## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.42 3.26
## 2 3.37 2.54
## 3 4.42 3.40
## 4 4.80 3.29
## 5 5.02 3.67
## 6 3.65 3.07
## 7 3.66 2.78
## 8 4.79 3.61
## 9 3.96 2.98
## 10 5.58 4.37menghitung rata-rata untuk 10 folds
mean_metric_polyC4 <- colMeans(metric_polyC4)
mean_metric_polyC4## rmse mae
## 4.366886 3.296632- Regresi polinomial derajat 5
set.seed(1501211)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_polyC5 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(horsepower,5,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_polyC5## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 5.82 4.18
## 2 3.89 2.91
## 3 4.37 3.41
## 4 3.34 2.60
## 5 3.80 2.97
## 6 4.08 3.00
## 7 4.09 2.91
## 8 5.10 3.70
## 9 4.29 3.34
## 10 5.02 3.85menghitung rata-rata 10 folds
mean_metric_polyC5 <- colMeans(metric_polyC5)
mean_metric_polyC5## rmse mae
## 4.380609 3.286917- Komparasi model regresi polinomial
## polynomZ rmse mae
## poly 1 4.875532 3.835955
## poly 2 4.334441 3.280627
## poly 3 4.286842 3.259572
## poly 4 4.366886 3.296632
## poly 5 4.380609 3.286917Berdasarkan hitungan dipoleh regresi polinomial berderajat 3 yang memiliki nilai rmse dan mae terkecil.
ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~poly(x,3,raw=T),
lty = 1, col = "red",se = F)+
theme_bw()
Regresi fungsi tangga (Visual)
- breaks = 2
model_tanggaC2 = lm(mpg ~ cut(horsepower,2),data=Auto)
summary(model_tanggaC2)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 2), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9879 -4.4012 -0.4012 3.6238 20.6121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.9879 0.3522 73.78 <2e-16 ***
## cut(horsepower, 2)(138,230] -11.5867 0.7520 -15.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.162 on 390 degrees of freedom
## Multiple R-squared: 0.3784, Adjusted R-squared: 0.3768
## F-statistic: 237.4 on 1 and 390 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,2),
lty = 1, col = "red",se = F)+
theme_bw()
- breaks = 3
model_tanggaC3 = lm(mpg ~ cut(horsepower,3),data=Auto)
summary(model_tanggaC3)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 3), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.1770 -3.6020 -0.2969 3.3730 19.4230
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.1770 0.3612 75.25 <2e-16 ***
## cut(horsepower, 3)(107,169] -9.8877 0.6752 -14.64 <2e-16 ***
## cut(horsepower, 3)(169,230] -13.8802 1.0854 -12.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.79 on 389 degrees of freedom
## Multiple R-squared: 0.4525, Adjusted R-squared: 0.4497
## F-statistic: 160.7 on 2 and 389 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,3),
lty = 1, col = "red",se = F)+
theme_bw()
- breaks = 4
model_tanggaC4 = lm(mpg ~ cut(horsepower,4),data=Auto)
summary(model_tanggaC4)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 4), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.0533 -2.7826 -0.2927 2.8593 17.5467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.0533 0.3581 81.12 <2e-16 ***
## cut(horsepower, 4)(92,138] -8.4506 0.5946 -14.21 <2e-16 ***
## cut(horsepower, 4)(138,184] -14.2707 0.7005 -20.37 <2e-16 ***
## cut(horsepower, 4)(184,230] -16.2004 1.2647 -12.81 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.001 on 388 degrees of freedom
## Multiple R-squared: 0.5926, Adjusted R-squared: 0.5894
## F-statistic: 188.1 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,4),
lty = 1, col = "red",se = F)+
theme_bw()
* breaks = 5
model_tanggaC5 = lm(mpg ~ cut(horsepower,5),data=Auto)
summary(model_tanggaC5)##
## 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-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "red",se = F)+
theme_bw()
- breaks = 6
model_tanggaC6 = lm(mpg ~ cut(horsepower,6),data=Auto)
summary(model_tanggaC6)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 6), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.0252 -2.9344 -0.1298 2.5172 14.5748
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.0252 0.4604 69.56 <2e-16 ***
## cut(horsepower, 6)(76.7,107] -8.0908 0.5947 -13.60 <2e-16 ***
## cut(horsepower, 6)(107,138] -12.2742 0.8109 -15.14 <2e-16 ***
## cut(horsepower, 6)(138,169] -16.9697 0.7850 -21.62 <2e-16 ***
## cut(horsepower, 6)(169,199] -18.3824 1.1187 -16.43 <2e-16 ***
## cut(horsepower, 6)(199,230] -19.3889 1.4821 -13.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.672 on 386 degrees of freedom
## Multiple R-squared: 0.6462, Adjusted R-squared: 0.6416
## F-statistic: 141 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,6),
lty = 1, col = "red",se = F)+
theme_bw()
Regresi fungsi tangga (Empiris)
set.seed(15012110)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
breaks <- 2:10
best_tanggaC <- map_dfr(breaks, function(i){
metric_tanggaC <- map_dfr(cross_val$splits,
function(x){
training <- Auto[x$in_id,]
training$horsepower <- cut(training$horsepower,i)
mod <- lm(mpg ~ horsepower,
data=training)
labs_x <- levels(mod$model[,2])
labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))
testing <- Auto[-x$in_id,]
horsepower_new <- cut(testing$horsepower,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
pred <- predict(mod,
newdata=list(horsepower=horsepower_new))
truth <- testing$mpg
data_eval <- na.omit(data.frame(truth,pred))
rmse <- mlr3measures::rmse(truth = data_eval$truth,
response = data_eval$pred
)
mae <- mlr3measures::mae(truth = data_eval$truth,
response = data_eval$pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_tanggaC
# menghitung rata-rata untuk 10 folds
mean_metric_tanggaC <- colMeans(metric_tanggaC)
mean_metric_tanggaC
})
best_tanggaC <- cbind(breaks=breaks,best_tanggaC)
# menampilkan hasil all breaks
best_tanggaC## breaks rmse mae
## 1 2 6.117025 4.785989
## 2 3 5.746927 4.500392
## 3 4 4.978185 3.776735
## 4 5 4.703610 3.567789
## 5 6 4.650068 3.533650
## 6 7 4.503009 3.360023
## 7 8 4.389183 3.352657
## 8 9 4.562461 3.448243
## 9 10 4.583278 3.436204diperoleh rata-rata 10 folds
mean_metric_tanggaC<-colMeans(best_tanggaC)
mean_metric_tanggaC<-mean_metric_tanggaC[-1]
mean_metric_tanggaC## rmse mae
## 4.914861 3.751298model dengan rmse terkecil
best_tanggaC %>% slice_min(rmse)## breaks rmse mae
## 1 8 4.389183 3.352657model dengan mae terkecil
best_tanggaC %>% slice_min(mae)## breaks rmse mae
## 1 8 4.389183 3.352657Diperoleh rmse terkecil pada regresi fungsi tangga dengan breaks 8
model_tanggaC8 = lm(mpg ~ cut(horsepower,8),data=Auto)
summary(model_tanggaC8)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 8), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.9471 -2.6757 -0.1533 2.4015 14.5529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9085 0.5771 58.76 <2e-16 ***
## cut(horsepower, 8)(69,92] -6.9614 0.6910 -10.07 <2e-16 ***
## cut(horsepower, 8)(92,115] -12.7551 0.7425 -17.18 <2e-16 ***
## cut(horsepower, 8)(115,138] -15.6656 1.1264 -13.91 <2e-16 ***
## cut(horsepower, 8)(138,161] -18.7799 0.8568 -21.92 <2e-16 ***
## cut(horsepower, 8)(161,184] -19.9735 1.1470 -17.41 <2e-16 ***
## cut(horsepower, 8)(184,207] -21.1228 1.7720 -11.92 <2e-16 ***
## cut(horsepower, 8)(207,230] -21.0085 1.5159 -13.86 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.433 on 384 degrees of freedom
## Multiple R-squared: 0.6832, Adjusted R-squared: 0.6774
## F-statistic: 118.3 on 7 and 384 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="grey") +
stat_smooth(method = "lm",
formula = y~cut(x,8),
lty = 1, col = "red",se = F)+
theme_bw()
Regresi Natural cubic spline (Visual)
- Regresi Natural Cubic Spline df = 3
attr(ns(Auto$horsepower, df=3),"knots")## 33.33333% 66.66667%
## 84 110model_splineh3 = lm(mpg ~ ns(horsepower, knots = c(84,110)),data=Auto)
summary(model_splineh3)##
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = c(84, 110)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.8308 -2.4068 -0.1672 2.2524 15.9184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.276 1.039 36.85 <2e-16 ***
## ns(horsepower, knots = c(84, 110))1 -21.311 1.054 -20.22 <2e-16 ***
## ns(horsepower, knots = c(84, 110))2 -35.274 2.313 -15.25 <2e-16 ***
## ns(horsepower, knots = c(84, 110))3 -19.669 1.435 -13.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.362 on 388 degrees of freedom
## Multiple R-squared: 0.6901, Adjusted R-squared: 0.6877
## F-statistic: 288 on 3 and 388 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="orange")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(84,110)),
lty = 1,se=F)
- Regresi Natural Cubic Spline df = 4
attr(ns(Auto$horsepower, df=4),"knots")## 25% 50% 75%
## 75.0 93.5 126.0model_splineh4 = lm(mpg ~ ns(horsepower, knots = c(75,93.5,126)),data=Auto)
summary(model_splineh4)##
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = c(75, 93.5, 126)),
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5781 -2.4975 -0.2753 2.0543 15.7782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.019 1.326 27.167 <2e-16
## ns(horsepower, knots = c(75, 93.5, 126))1 -15.110 1.308 -11.554 <2e-16
## ns(horsepower, knots = c(75, 93.5, 126))2 -20.243 1.290 -15.692 <2e-16
## ns(horsepower, knots = c(75, 93.5, 126))3 -26.848 3.019 -8.892 <2e-16
## ns(horsepower, knots = c(75, 93.5, 126))4 -21.761 1.569 -13.874 <2e-16
##
## (Intercept) ***
## ns(horsepower, knots = c(75, 93.5, 126))1 ***
## ns(horsepower, knots = c(75, 93.5, 126))2 ***
## ns(horsepower, knots = c(75, 93.5, 126))3 ***
## ns(horsepower, knots = c(75, 93.5, 126))4 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.328 on 387 degrees of freedom
## Multiple R-squared: 0.6957, Adjusted R-squared: 0.6925
## F-statistic: 221.2 on 4 and 387 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="orange")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(75,93.5,126)),
lty = 1,se=F)
- Regresi Natural Cubic Spline df = 5
attr(ns(Auto$horsepower, df=5),"knots")## 20% 40% 60% 80%
## 72 88 100 140model_splineh5 = lm(mpg ~ ns(horsepower, knots = c(72,88,100,140)),data=Auto)
summary(model_splineh5)##
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = c(72, 88, 100, 140)),
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.9340 -2.5874 -0.1347 2.1908 15.5894
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 34.841 1.442 24.165
## ns(horsepower, knots = c(72, 88, 100, 140))1 -10.739 1.468 -7.314
## ns(horsepower, knots = c(72, 88, 100, 140))2 -14.381 1.794 -8.015
## ns(horsepower, knots = c(72, 88, 100, 140))3 -21.688 1.437 -15.094
## ns(horsepower, knots = c(72, 88, 100, 140))4 -23.028 3.299 -6.979
## ns(horsepower, knots = c(72, 88, 100, 140))5 -21.381 1.626 -13.149
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## ns(horsepower, knots = c(72, 88, 100, 140))1 1.51e-12 ***
## ns(horsepower, knots = c(72, 88, 100, 140))2 1.33e-14 ***
## ns(horsepower, knots = c(72, 88, 100, 140))3 < 2e-16 ***
## ns(horsepower, knots = c(72, 88, 100, 140))4 1.30e-11 ***
## ns(horsepower, knots = c(72, 88, 100, 140))5 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.311 on 386 degrees of freedom
## Multiple R-squared: 0.6989, Adjusted R-squared: 0.695
## F-statistic: 179.2 on 5 and 386 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="orange")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(72,88,100,140)),
lty = 1,se=F)
- Regresi Natural Cubic Spline df = 6
attr(ns(Auto$horsepower, df=6),"knots")## 16.66667% 33.33333% 50% 66.66667% 83.33333%
## 70.0 84.0 93.5 110.0 150.0model_splineh6 = lm(mpg ~ ns(horsepower, knots = c(70,84,93.5,110,150)),data=Auto)
summary(model_splineh6)##
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = c(70, 84, 93.5, 110,
## 150)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.9491 -2.6183 -0.1595 2.3508 15.1349
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 34.738 1.509 23.021
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))1 -8.210 1.594 -5.149
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))2 -13.046 1.835 -7.108
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))3 -14.577 1.886 -7.730
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))4 -22.802 1.624 -14.039
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))5 -22.758 3.512 -6.480
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))6 -20.849 1.742 -11.967
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))1 4.18e-07 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))2 5.76e-12 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))3 9.50e-14 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))4 < 2e-16 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))5 2.81e-10 ***
## ns(horsepower, knots = c(70, 84, 93.5, 110, 150))6 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.302 on 385 degrees of freedom
## Multiple R-squared: 0.7009, Adjusted R-squared: 0.6962
## F-statistic: 150.4 on 6 and 385 DF, p-value: < 2.2e-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="orange")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(70,84,93.5,110,150)),
lty = 1,se=F)
Regresi Natural Cubic Spline (Empiris)
library(splines)
set.seed(1501211044)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 2:10
best_splineC <- map_dfr(df, function(i){
metric_splineC <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(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_splineC
# menghitung rata-rata untuk 10 folds
mean_metric_splineC <- colMeans(metric_splineC)
mean_metric_splineC
}
)menampilkan hasil all breaks
best_splineC <- cbind(df=df,best_splineC)
best_splineC## df rmse mae
## 1 2 4.325352 3.265504
## 2 3 4.346849 3.279236
## 3 4 4.318819 3.261885
## 4 5 4.315448 3.261790
## 5 6 4.314716 3.252722
## 6 7 4.288756 3.223193
## 7 8 4.286997 3.237556
## 8 9 4.287206 3.222053
## 9 10 4.283672 3.211760rata-rata untuk 10 folds
mean_metric_C<-colMeans(best_splineC[-1])
mean_metric_C## rmse mae
## 4.307535 3.246189model terbaik berdasarkan rmse terkecil
best_splineC %>% slice_min(rmse)## df rmse mae
## 1 10 4.283672 3.21176model terbaik berdasarkan mae terkecil
best_splineC %>% slice_min(mae)## df rmse mae
## 1 10 4.283672 3.21176Berdasarkan hitungan dipilih model natural cubic spline dengan df = 10 yang memiliki niai rmse dan mae terkecil. Sehingga digunakan df=10 untuk regresi natural cubic spline terbaik, untuk knots berdasarkan df=10 dihitung menggunakan komputer.
attr(ns(Auto$horsepower, df=10),"knots")## 10% 20% 30% 40% 50% 60% 70% 80% 90%
## 67.0 72.0 80.0 88.0 93.5 100.0 110.0 140.0 157.7model_splineC10 = lm(mpg ~ ns(horsepower, knots = c(67,72,80,88,93.5,100,110,140,157.7)),data=Auto)
summary(model_splineC10)##
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = c(67, 72, 80, 88, 93.5,
## 100, 110, 140, 157.7)), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5003 -2.4956 -0.2086 2.2577 14.6400
##
## Coefficients:
## Estimate
## (Intercept) 32.462
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))1 -4.600
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))2 -3.007
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))3 -8.663
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))4 -7.170
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))5 -14.722
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))6 -8.446
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))7 -16.947
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))8 -21.715
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))9 -14.384
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))10 -22.326
## Std. Error
## (Intercept) 1.680
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))1 1.789
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))2 2.555
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))3 2.043
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))4 2.188
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))5 2.125
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))6 2.467
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))7 2.080
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))8 1.891
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))9 4.116
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))10 1.946
## t value
## (Intercept) 19.324
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))1 -2.572
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))2 -1.177
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))3 -4.239
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))4 -3.278
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))5 -6.927
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))6 -3.424
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))7 -8.148
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))8 -11.483
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))9 -3.494
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))10 -11.476
## Pr(>|t|)
## (Intercept) < 2e-16
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))1 0.010495
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))2 0.239996
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))3 2.82e-05
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))4 0.001142
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))5 1.84e-11
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))6 0.000684
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))7 5.38e-15
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))8 < 2e-16
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))9 0.000531
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))10 < 2e-16
##
## (Intercept) ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))1 *
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))2
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))3 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))4 **
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))5 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))6 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))7 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))8 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))9 ***
## ns(horsepower, knots = c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7))10 ***
## ---
## 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-16ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="orange")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = c(67,72,80,88,93.5,100,110,140,157.7)),
lty = 1,se=F)
- Regresi natural cubic spline 9 knots
set.seed(2021)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_splinek9 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(67,72,80,88,93.5,100,110,140,157.7)),
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_splinek9## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.69 2.95
## 2 3.90 3.03
## 3 4.46 3.23
## 4 4.40 3.39
## 5 5.53 4.02
## 6 3.52 2.56
## 7 5.04 3.55
## 8 4.54 3.39
## 9 4.37 3.20
## 10 3.61 2.63rata-rata untuk 10 fold
mean_metric_splinek9 <- colMeans(metric_splinek9)
mean_metric_splinek9## rmse mae
## 4.305265 3.195052Komparasi model
## model rmse mae
## Linear 4.883887 3.835224
## Poly 1 4.875532 3.835955
## Poly 2 4.334441 3.280627
## Poly 3 4.286842 3.259572
## Poly 4 4.366886 3.296632
## Poly5 4.380609 3.286917
## Tangga 4.914861 3.751298
## Spline df=10 4.307535 3.246189
## Spline 9 knots 4.305265 3.195052z %>% slice_min(rmse)## model rmse mae
## mean_metric_polyC3 Poly 3 4.286842 3.259572z %>% slice_min(mae)## model rmse mae
## mean_metric_splinek9 Spline 9 knots 4.305265 3.195052Berdasarkan hasil tersebut, diperoleh model regresi polynomial berderajat 3 yang memiliki rmse terkecil dan regresi natural cubic spline dengan 9 knots yang memiliki nilai mae terkecil 