folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()
set.seed(888)
for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(horsepower,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

poly_model <- lm(mpg ~ poly(horsepower, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$Horsepower^', row_number()-1,'$'))) %>%
  kable(digits = 3)

poly6 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('Horsepower') + ylab('Miles Per Gallon') +
  geom_point(aes(horsepower, mpg, col = 'blue')) +
  geom_line(aes(horsepower, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
poly6

folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()
for (k in 1:10) {
  for (i in 1:10) {
    Auto$cut <- cut(Auto$horsepower, i+1)
    model <- lm(mpg ~ horsepower*cut, data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
    Auto$cut <- NULL
  }
}
errors <- apply(errors, 2, mean)

data_frame(RMSE = errors) %>%
  mutate(Steps = row_number() +1) %>%
  ggplot(aes(Steps, RMSE, fill = Steps == (which.min(errors) + 1))) + 
  geom_col() +
  scale_x_continuous(breaks = 2:(length(errors)+1)) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

#dari gambar dapat dilihat bahwa untuk horsepower diperlukan 7 potongan 

min_error <- which.min(errors) + 1
step_model <- lm(mpg ~ horsepower *cut(horsepower, min_error), data = Auto)
tidy(step_model) %>%
  kable(digits = 3)

pc7 <- Auto %>%
  mutate(Predictions = predict(step_model, Auto)) %>%
  ggplot() + xlab('Horse Power') + ylab('Miled Per Gallon') + 
  geom_point(aes(horsepower, mpg, col = 'blue')) +
  geom_line(aes(horsepower, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))

pc7

set.seed(888)
spline_model <- train(mpg ~ horsepower, data
                      = Auto,
                      method = 'gamSpline',
                      trControl = trainControl(method = 'cv', number = 10),
                      tuneGrid = expand.grid(df = seq(1, 12, 1)))

disp_lm <- train(mpg ~ horsepower, data = Auto, method = 'lm',
                 trControl = trainControl(method = 'cv', number = 10))
plot(spline_model)

anova(disp_lm$finalModel, 
      spline_model$finalModel) %>%
  tidy %>%
  mutate(model = c('Linear', 'Spline'),
         percent = sumsq/lag(rss) * 100) %>% 
  select(model, rss, percent, sumsq, p.value) %>%
  kable(digits = 3)

splineK12 <- Auto %>%
  mutate(pred = predict(spline_model, Auto)) %>%
  ggplot() + 
  geom_point(aes(horsepower, mpg, col = 'blue')) +
  geom_line(aes(horsepower, pred, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00')) +
  labs(x = 'Horse Power', y = 'Miles Per Gallon', title = 'Spline Relationship') +
  theme(legend.position = 'top')

splineK12

folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()
set.seed (888)
for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(acceleration,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)
errors

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

poly_model <- lm(mpg ~ poly(acceleration, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$acceleration^', row_number()-1,'$'))) %>%
  kable(digits = 3)

polys1 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('acceleration') + ylab('Miles Per Gallon') +
  geom_point(aes(acceleration, mpg, col = 'blue')) +
  geom_line(aes(acceleration, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
polys1

folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()
set.seed (888)
for (k in 1:10) {
  for (i in 1:10) {
    Auto$cut <- cut(Auto$acceleration, i+1)
    model <- lm(mpg ~ acceleration*cut, data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
    Auto$cut <- NULL
  }
}
errors<- apply(errors, 2, mean)
errors

data_frame(RMSE = errors) %>%
  mutate(Steps = row_number() +1) %>%
  ggplot(aes(Steps, RMSE, fill = Steps == (which.min(errors) + 1))) + 
  geom_col() +
  scale_x_continuous(breaks = 2:(length(errors)+1)) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

#dari gambar dapat dilihat bahwa untuk acceleration diperlukan 10 potongan 

min_error <- which.min(errors) + 1
step_model <- lm(mpg ~ acceleration *cut(acceleration, min_error), data = Auto)
tidy(step_model) %>%
  kable(digits = 3)

pc1 <- Auto %>%
  mutate(Predictions = predict(step_model, Auto)) %>%
  ggplot() + xlab('acceleration') + ylab('Miled Per Gallon') + 
  geom_point(aes(acceleration, mpg, col = 'blue')) +
  geom_line(aes(acceleration, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
pc1

spline_model <- train(mpg ~ acceleration, data
                      = Auto,
                      method = 'gamSpline',
                      trControl = trainControl(method = 'cv', number = 10),
                      tuneGrid = expand.grid(df = seq(1, 12, 1)))
disp_lm <- train(mpg ~ acceleration, data = Auto, method = 'lm',
                 trControl = trainControl(method = 'cv', number = 10))
disp_lm

plot(spline_model)

anova(disp_lm$finalModel, 
      spline_model$finalModel) %>%
  tidy %>%
  mutate(model = c('Linear', 'Spline'),
         percent = sumsq/lag(rss) * 100) %>% 
  select(model, rss, percent, sumsq, p.value) %>%
  kable(digits = 3)

splineK4 <- Auto %>%
  mutate(pred = predict(spline_model, Auto)) %>%
  ggplot() + 
  geom_point(aes(acceleration, mpg, col = 'blue')) +
  geom_line(aes(acceleration, pred, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00')) +
  labs(x = 'acceleration', y = 'Miles Per Gallon', title = 'Spline Relationship') +
  theme(legend.position = 'top')

splineK4

#regresi polinomial dengan CV-10
folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()

for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(displacement,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

#model terbaik polinomial adalah derajat 5 dari grafik tersebut

poly_model <- lm(mpg ~ poly(displacement, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$Horsepower^', row_number()-1,'$'))) %>%
  kable(digits = 3)

#namun ketika melihat nilai p.value nya yang signifikan hanya polinomial derajat 1,2,
#berarti menunjukkan bahwa antara mpg dan weight memiliki hubungan non-linier.

poly5 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('displacement') + ylab('Miles Per Gallon') +
  geom_point(aes(displacement, mpg, col = 'blue')) +
  geom_line(aes(displacement, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
poly5

set.seed(888)
folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()
for (k in 1:10) {
  for (i in 1:10) {
    Auto$cut <- cut(Auto$displacement, i+1)
    model <- lm(mpg ~ displacement*cut, data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
    Auto$cut <- NULL
  }
}
errors <- apply(errors, 2, mean)

data_frame(RMSE = errors) %>%
  mutate(Steps = row_number() +1) %>%
  ggplot(aes(Steps, RMSE, fill = Steps == (which.min(errors) + 1))) + 
  geom_col() +
  scale_x_continuous(breaks = 2:(length(errors)+1)) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

#dari gambar dapat dilihat bahwa untuk horsepower diperlukan 7 potongan 

min_error <- which.min(errors) + 1
step_model <- lm(mpg ~ displacement *cut(displacement, min_error), data = Auto)
tidy(step_model) %>%
  kable(digits = 3)

pc5 <- Auto %>%
  mutate(Predictions = predict(step_model, Auto)) %>%
  ggplot() + xlab('displacement') + ylab('Miled Per Gallon') + 
  geom_point(aes(displacement, mpg, col = 'blue')) +
  geom_line(aes(displacement, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
pc5

spline_model <- train(mpg ~ displacement, data
                      = Auto,
                      method = 'gamSpline',
                      trControl = trainControl(method = 'cv', number = 10),
                      tuneGrid = expand.grid(df = seq(1, 12, 1)))
disp_lm <- train(mpg ~ displacement, data = Auto, method = 'lm',
                 trControl = trainControl(method = 'cv', number = 10))
plot(spline_model)

anova(disp_lm$finalModel, 
      spline_model$finalModel) %>%
  tidy %>%
  mutate(model = c('Linear', 'Spline'),
         percent = sumsq/lag(rss) * 100) %>% 
  select(model, rss, percent, sumsq, p.value) %>%
  kable(digits = 3)

splineK12 <- Auto %>%
  mutate(pred = predict(spline_model, Auto)) %>%
  ggplot() + 
  geom_point(aes(displacement, mpg, col = 'blue')) +
  geom_line(aes(displacement, pred, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00')) +
  labs(x = 'displacement', y = 'Miles Per Gallon', title = 'Spline Relationship') +
  theme(legend.position = 'top')
splineK12

set.seed(888)
folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()

for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(horsepower,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)
errors

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

poly_model <- lm(mpg ~ poly(horsepower, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$Horsepower^', row_number()-1,'$'))) %>%
  kable(digits = 3)

poly6 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('Horsepower') + ylab('Miles Per Gallon') +
  geom_point(aes(horsepower, mpg, col = 'blue')) +
  geom_line(aes(horsepower, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
poly6

cv.poly3 <- function(dataset){
  set.seed(888)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly3 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,3,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-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)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
}

cv.poly3(Auto)

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

#perbandingan antara polinomial derajat 3 dan 6
ggplot(Auto,aes(x=horsepower,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="lm",
              formula=y~poly(x,6,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="blue",se=F)+
    stat_smooth(method="lm",
              formula=y~poly(x,3,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="red",se=F)+
  scale_color_manual(values = c("Regresi Polinomial"="blue","Regresi Polinomial"="red"))

set.seed(888)

cv.pc <- function(dataset){
  set.seed(888)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 3:12

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[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 <- dataset[-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)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
basedonrmse <- as.data.frame(best_tangga %>% slice_min(rmse))
basedonmae <- as.data.frame(best_tangga %>% slice_min(mae))
return(rbind(basedonrmse,basedonmae))
}

cv.pc(Auto)

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

pc12 <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,12), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

plot_grid(pc8, pc12, labels = c("pc8", "pc12"), label_size = 6)

set.seed(888)
cross.val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 1:12
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- 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.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-data.frame(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                 best.spline3 %>% slice_min(mae))) #berdasarkan mae

basis

df_7<-attr(ns(Auto$horsepower, df=7),"knots")
df_7

df_10<-attr(ns(Auto$horsepower, df=10),"knots")
df_10

set.seed(888)

cross.val <- vfold_cv(Auto,v=10,strata = "mpg")

metric.spline3.ns6 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(68.85714, 79.71429, 89.57143, 98.00000, 113.57143,150.00000)),
         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)
}
)
mean.metric.spline3.ns6 <- colMeans(metric.spline3.ns6)# menghitung rata-rata untuk 10 folds



metric.spline3.ns9 <- 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)
}
)
mean.metric.spline3.ns9 <- colMeans(metric.spline3.ns9)# menghitung rata-rata untuk 10 folds



nilai.cv.ns <- rbind("NCS6" = mean.metric.spline3.ns6,
                     "NCS9" = mean.metric.spline3.ns9)

nama.model.ncs <- c("6 - NCS (df=7)",
                    "9 - NCS (df=10)")

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs

nc6 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="darkgrey") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(68.85714, 79.71429, 89.57143, 98.00000, 113.57143,150.00000)), col = "blue", se = F) +
  theme_bw()


nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="darkgrey") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(67,72,80,88,93.5,100,110,140,157.7)), col = "blue", se = F)+
  theme_bw()

plot_grid(nc6, nc9, labels = c("nc6", "nc9"), label_size = 6)

perbandingan.model <- rbind(cv.poly3(Auto), cv.pc(Auto)[1,-1], model.ncs[2,-1])
perbandingan.model$metode <- c("Polinomial Ordo 3","Piecewise Constant 12","NCS 9 Knot")
perbandingan.model

plot_grid(poly3, pc12, nc9, labels = c("poly3","pc12", "nc9"), label_size = 6)

ggplot(Auto,aes(x=horsepower,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="lm",
              formula=y~poly(x,3,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="blue",se=F)+
  stat_smooth(method="lm",
              formula=y~cut(x,12),
            lty=1,aes(col="Piecewise Constant"),col="red",se=F)+
   stat_smooth(method="lm",
              formula=y~ns(x,df=9),
            lty=1,aes(col="Natural Cubic Spline"),col="black",se=F)+
  scale_color_manual(values = c("Regresi Polinomial"="blue","Piecewise Constant"="red","Natural Cubic Spline"="black"))

set.seed(888)
folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()

for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(acceleration,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)
errors

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

poly_model <- lm(mpg ~ poly(acceleration, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$acceleration^', row_number()-1,'$'))) %>%
  kable(digits = 3)

poly1 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('acceleration') + ylab('Miles Per Gallon') +
  geom_point(aes(acceleration, mpg, col = 'blue')) +
  geom_line(aes(acceleration, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
poly1

set.seed(888)

cv.pc <- function(dataset){
  set.seed(888)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 3:12

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[x$in_id,]
training$acceleration <- cut(training$acceleration,i)

mod <- lm(mpg ~ acceleration,
         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 <- dataset[-x$in_id,]

acceleration_new <- cut(testing$acceleration,c(labs_x_breaks[1,1],labs_x_breaks[,2]))

pred <- predict(mod,
               newdata=list(acceleration=acceleration_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)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
basedonrmse <- as.data.frame(best_tangga %>% slice_min(rmse))
basedonmae <- as.data.frame(best_tangga %>% slice_min(mae))
return(rbind(basedonrmse,basedonmae))
}

cv.pc(Auto)

pc3 <- ggplot(Auto,aes(x=acceleration, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,3), 
               lty = 1, col = "red",se = F)+
  theme_bw()

pc3

set.seed(888)
cross.val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 1:12
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- map_dfr(cross.val$splits,
                  function(x){
                  mod <- lm(mpg ~ ns(acceleration,df=i),data=Auto[x$in_id,])
                  pred <- predict(mod,newdata=Auto[-x$in_id,])
                  truth <- Auto[-x$in_id,]$mpg
                  rmse <- mlr3measures::rmse(truth = truth,response = pred)
                  mae <- mlr3measures::mae(truth = truth,response = pred)
                  metric <- c(rmse,mae)
                  names(metric) <- c("rmse","mae")
                  return(metric)
                }
      )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-data.frame(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                 best.spline3 %>% slice_min(mae))) #berdasarkan mae

basis

df_5<-attr(ns(Auto$acceleration, df=5),"knots")
df_5

set.seed(888)

cross.val <- vfold_cv(Auto,v=10,strata = "mpg")

metric.spline3.ns5 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(acceleration, knots = c(13.42, 14.80, 16.00, 17.70 )),
         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)
}
)
mean.metric.spline3.ns5 <- colMeans(metric.spline3.ns5)# menghitung rata-rata untuk 10 folds



nilai.cv.ns <- mean.metric.spline3.ns5
                

nama.model.ncs <- "4 - NCS (df=5)"
                    

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs

nc4 <- ggplot(Auto, aes(x = acceleration, y = mpg)) + 
  geom_point(alpha = 0.55, color="darkgrey") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(13.42, 14.80, 16.00, 17.70 )), col = "blue", se = F)+
  theme_bw()

nc4

perbandingan.model <- rbind(5.832256, cv.pc(Auto)[1,-1], model.ncs[2,-1])
perbandingan.model$metode <- c("Polinomial Ordo 1","Piecewise Constant 3","NCS 4
                               Knot")
perbandingan.model

plot_grid(poly1, pc3, nc4, labels = c("poly1","pc3", "nc4"), label_size = 6)

ggplot(Auto,aes(x=acceleration,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="lm",
              formula=y~poly(x,1,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="blue",se=F)+
  stat_smooth(method="lm",
              formula=y~cut(x,3),
            lty=1,aes(col="Piecewise Constant"),col="red",se=F)+
   stat_smooth(method="lm",
              formula=y~ns(x,df=5),
            lty=1,aes(col="Natural Cubic Spline"),col="black",se=F)+
  scale_color_manual(values = c("Regresi Polinomial"="blue","Piecewise Constant"="red","Natural Cubic Spline"="black"))

set.seed(888)
folds <- sample(rep(1:10, nrow(Auto)/10))
errors <- matrix(NA, 10, 10)
models <- list()

for (k in 1:10) {
  for (i in 1:10) {
    model <- lm(mpg ~ poly(displacement,i), data = Auto[folds != i,])
    pred <- predict(model, Auto[folds == i,])
    resid <- (Auto$mpg[folds==i] - pred)^2
    errors[k, i] <- sqrt(mean(resid))
  }
}
errors <- apply(errors, 2, mean)
errors

data_frame(RMSE = errors) %>%
  mutate(Polynomial = row_number()) %>%
  ggplot(aes(Polynomial, RMSE, fill = Polynomial == which.min(errors))) + 
  geom_col() +
  scale_x_continuous(breaks = 1:10) +
  guides(fill = FALSE) +
  coord_cartesian(ylim = c(min(errors), max(errors)))

poly_model <- lm(mpg ~ poly(displacement, which.min(errors)), data = Auto)
tidy(poly_model) %>%
  mutate(term = ifelse(row_number() == 1, 'Intercept',
                       paste0('$displacement^', row_number()-1,'$'))) %>%
  kable(digits = 3)

polys6 <- Auto %>%
  mutate(Predictions = predict(poly_model, Auto)) %>%
  ggplot() + xlab('displacement') + ylab('Miles Per Gallon') +
  geom_point(aes(displacement, mpg, col = 'blue')) +
  geom_line(aes(displacement, Predictions, col = 'goldenrod2'), size = 1.5) +
  scale_color_manual(name = 'Value Type',
                     labels = c('Observed', 'Predicted'),
                     values = c('#56B4E9', '#E69F00'))
polys6

cv.poly4 <- function(dataset){
  set.seed(888)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly4 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(displacement,4,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-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)
      }
  )
  poly4<- colMeans(metric_poly4) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly4)
}

cv.poly4(Auto)

poly4 <- ggplot(Auto,aes(x=displacement, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,4,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()
poly4

#perbandingan antara polinomial derajat 4 dan 6
ggplot(Auto,aes(x=displacement,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="lm",
              formula=y~poly(x,6,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="blue",se=F)+
    stat_smooth(method="lm",
              formula=y~poly(x,4,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="red",se=F)+
  scale_color_manual(values = c("Regresi Polinomial"="blue","Regresi Polinomial"="red"))

set.seed(888)

cv.pc <- function(dataset){
  set.seed(888)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 3:12

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[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 <- dataset[-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)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
basedonrmse <- as.data.frame(best_tangga %>% slice_min(rmse))
basedonmae <- as.data.frame(best_tangga %>% slice_min(mae))
return(rbind(basedonrmse,basedonmae))
}

cv.pc(Auto)

pc9 <- 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 = "red",se = F)+
  theme_bw()

pc9

set.seed(888)
cross.val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 1:12
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- 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.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-data.frame(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                 best.spline3 %>% slice_min(mae))) #berdasarkan mae

basis

df_12<-attr(ns(Auto$displacement, df=12),"knots")
df_12

set.seed(888)

cross.val <- vfold_cv(Auto,v=10,strata = "mpg")

metric.spline3.ns12 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(acceleration, knots = c(12.00000, 13.00000, 13.77500, 14.50000, 14.90000, 15.50000, 16.00000, 16.50000, 17.02500 )),
         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)
}
)
mean.metric.spline3.ns12 <- colMeans(metric.spline3.ns12)# menghitung rata-rata untuk 10 folds



nilai.cv.ns <- mean.metric.spline3.ns12
                

nama.model.ncs <- "11 - NCS (df=12)"
                    

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs

nc11 <- ggplot(Auto, aes(x = displacement, y = mpg)) + 
  geom_point(alpha = 0.55, color="darkgrey") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(12.00000, 13.00000, 13.77500, 14.50000, 14.90000, 15.50000, 16.00000, 16.50000, 17.02500 )), col = "blue", se = F)+
  theme_bw()

nc11

perbandingan.model <- rbind(cv.poly4(Auto), cv.pc(Auto)[1,-1], model.ncs[2,-1])
perbandingan.model$metode <- c("Polinomial Ordo 4","Piecewise Constant 9","NCS 11 Knot")
perbandingan.model

plot_grid(poly4, pc9, nc11, labels = c("poly4","pc9", "nc11"), label_size = 6)

ggplot(Auto,aes(x=displacement,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="lm",
              formula=y~poly(x,4,raw=T),
            lty=1,aes(col="Regresi Polinomial"),col="blue",se=F)+
  stat_smooth(method="lm",
              formula=y~cut(x,9),
            lty=1,aes(col="Piecewise Constant"),col="red",se=F)+
   stat_smooth(method = "lm", formula = y~ns(x, knots = c(12.00000, 13.00000, 13.77500, 14.50000, 14.90000, 15.50000, 16.00000, 16.50000, 17.02500 )), col = "black", se = F)+
  scale_color_manual(values = c("Regresi Polinomial"="blue","Piecewise Constant"="red","Natural Cubic Spline"="black"))

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

pred_data <- broom::augment(model_sms)

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

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

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

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

pred_data <- broom::augment(model_sms)

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

model_sms2 <- with(data = Auto,smooth.spline(acceleration,mpg))
model_sms2                  

pred_data2 <- broom::augment(model_sms2)

smooth2 <- ggplot(pred_data2,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("acceleration")+
  ylab("mpg")+
  theme_bw()

smooth2

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

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

model_sms21 <- with(data = Auto,smooth.spline(acceleration,mpg,df=3))
model_sms21                  

pred_data21 <- broom::augment(model_sms21)

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

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

pred_data3 <- broom::augment(model_sms)

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

smooth3

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

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

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

pred_data31 <- broom::augment(model_sms31)

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

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

model_loess_span <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~horsepower,
                     data = Auto,span=.$span)))

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

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

set.seed(888)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")

span <- seq(0.1,0.5,length.out=50)

best_loess1 <- map_dfr(span, function(i){
metric_loess <- map_dfr(cross_val$splits,
    function(x){
mod <- loess(mpg ~ horsepower,span = i,
         data=Auto[x$in_id,])
pred <- predict(mod,
               newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg

data_eval <- na.omit(data.frame(pred=pred,
                                truth=truth))


rmse <- mlr3measures::rmse(truth = data_eval$truth,
                           response = data_eval$pred
                           )
mae <- mlr3measures::mae(truth = data_eval$truth,
                           response = data_eval$pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)

head(metric_loess,20)

# menghitung rata-rata untuk 10 folds
mean_metric_loess <- colMeans(metric_loess)

mean_metric_loess
}
)

best_loess1 <- cbind(span=span,best_loess1)
# menampilkan hasil all breaks
best_loess1

#berdasarkan rmse
best_loess1 %>% slice_min(rmse)

#berdasarkan mae
best_loess1 %>% slice_min(mae)

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

ggplot(Auto,aes(x=horsepower,y=mpg))+
  geom_point(alpha=0.55,color="darkgrey")+
  stat_smooth(method="loess",
              formula=y~x,span = 0.1081633,col="blue",lty=1,se=F)+
  stat_smooth(method="loess",
              formula=y~x,span = 0.75,col="red",lty=1,se=F)+
  scale_color_manual(values = c("LOESS 0.1081633"="blue","OESS 0.75"="red"))

model_loess <- loess(mpg~acceleration,
                     data = Auto)
summary(model_loess)

model_loess_span <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~acceleration,
                     data = Auto,span=.$span)))

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

loess2 <- ggplot(Auto, aes(acceleration,mpg)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.75,
              col="blue",
              lty=1,
              se=F)
loess2

set.seed(888)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")

span <- seq(0.1,0.5,length.out=50)

best_loess2 <- map_dfr(span, function(i){
metric_loess <- map_dfr(cross_val$splits,
    function(x){
mod <- loess(mpg ~ acceleration,span = i,
         data=Auto[x$in_id,])
pred <- predict(mod,
               newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg

data_eval <- na.omit(data.frame(pred=pred,
                                truth=truth))


rmse <- mlr3measures::rmse(truth = data_eval$truth,
                           response = data_eval$pred
                           )
mae <- mlr3measures::mae(truth = data_eval$truth,
                           response = data_eval$pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)

head(metric_loess,20)

# menghitung rata-rata untuk 10 folds
mean_metric_loess <- colMeans(metric_loess)

mean_metric_loess
}
)

best_loess2 <- cbind(span=span,best_loess2)
# menampilkan hasil all breaks
best_loess2

#berdasarkan rmse
best_loess2 %>% slice_min(rmse)

#berdasarkan mae
best_loess2 %>% slice_min(mae)

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