Non Linear Model
Soal
Data yang digunakan diperoleh dari package ISLR. Gunakan model-model non-linear pada data Auto dengan respon mpg dan pilih tiga kolom untuk dijadikan peubah prediktor.
Apakah ada bukti untuk hubungan non-linear terhadap peubah-peubah yang dipilih? Buatlah visualisasi untuk mendukung jawaban Anda!
Tentukan model non-linear terbaik untuk masing-masing pasangan peubah (secara visual dan/atau secara empiris)
Menjalankan Library
Library yang digunakan pada pembahasan soal kali ini antara lain:
library(ISLR)
library(DT)
library(tidyverse)
library(splines)
library(rsample)Data yang digunakan pada pembahasan kali ini adalah data Auto dari package ISLR yang terdiri dari 392 observasi dengan 9 peubah. Namun, pada pembahasan soal kali ini, peubah yang digunakan ada 4, yaitu:
- mpg (mile per gallon)
- horsepower
- weight
- displacement
Pemodelan horsepower dan mpg
Menampilkan Data yang Digunakan
attach(Auto)
datatable(cbind(mpg, horsepower))plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Scatter Plot horsepower vs mpg")
Grafik di atas merupakan scatter plot dari 392 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah keuatan mesin kendaraan (horsepower), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).
Secara grafik, terlihat bahwa pola data menunjukkan ke pola nonlinear. Untuk menemukan solusi optimal dari data kita bisa menggunakan Model Regresi Polinomial, Fungsi Tangga (Piecewise Constant) dan Natural Cubic Splines.
Regresi Polinomial
Secara Empiris
Cross Validation Data Auto menggunakan Model Regresi Polinomial
- Regresi Polinomial Derajat 2
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(horsepower,2,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly2)- Regresi Polinomial Derajat 3
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly3 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(horsepower,3,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly3)- Regresi Polinomial Derajat 4
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly4 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(horsepower,4,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly4)Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.
poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)
cvpoly <- rbind (poly2,poly3,poly4)
cvpoly## rmse mae
## poly2 4.382221 3.287275
## poly3 4.403984 3.300528
## poly4 4.413418 3.316362Berdasarkan hasil rata-rata rsme dan mae regresi polinomial diatas, terlihat bahwa Regresi Polinomial Derajat 2 memiliki nilai rsme dan mae paling kecil. Artinya, Regresi Polinomial Derajat 2 lebih baik dibandingkan Regresi Polinomial Derajat 3 dan 4.
# menampilkan model
nilai.cvpoly <- rbind("Poly2" = poly2)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 2")
model_poly_h <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_h## Model rmse mae
## Poly2 Regresi Polinomial Derajat 2 4.382221 3.287275Secara Visualisasi
par(mfrow=c(2,2))
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Scatter Plot horsepower vs mpg")
polinom_d2 <- lm(mpg ~ poly(horsepower,2))
prediksi <- predict(polinom_d2, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 2")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d3 <- lm(mpg ~ poly(horsepower,3))
prediksi <- predict(polinom_d3, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 3")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d4 <- lm(mpg ~ poly(horsepower,4))
prediksi <- predict(polinom_d4, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 4")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
Fungsi Tangga (Piecewise Constant)
Secara Empiris
Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
breaks <- 3:12
set.seed(1501211024)
best_tangga <- map_dfr(breaks, function(i){
metric_tangga <- 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_tangga
# menghitung rata-rata untuk 8 folds
mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)
# menampilkan semua jumlah interval
best_tangga## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 5.74 4.52
## 2 5.02 3.82
## 3 4.73 3.60
## 4 4.62 3.53
## 5 4.55 3.39
## 6 4.49 3.43
## 7 4.57 3.50
## 8 4.56 3.44
## 9 4.50 3.35
## 10 4.50 3.36# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)## breaks rmse mae
## 1 8 4.487083 3.431063
## 2 11 4.504833 3.353963Berdasarkan rmse jumlah interval terbaik adalah 8, dan berdasarkan mae jumlah interval terbaik adalah 11. Karena selisih mae pada breaks 8 dan 11 kecil, maka mae dianggap sama, sehingg pada kasus ini breaks 8 atau Fungsi Tangga dengan 7 knots adalah interval terbaik.
Secara Visualiasi
par(mfrow=c(1,2))
tangga8 <- lm(mpg ~ cut(horsepower,7))
prediksi <- predict(tangga8, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 7 knots")
ix <- sort(horsepower, index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l",lwd=2)
tangga11 <- lm(mpg ~ cut(horsepower,10))
prediksi <- predict(tangga11, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 10 knots")
ix <- sort(horsepower, index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l",lwd=2)
Berdasrkan grafik, juga terlihat bahwa Model Fungsi Tangga dengan 7 knots lebih baik dari 10 knots.
# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 7 knots")
model_tangga_h <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_h ## Model rmse mae
## tangga Fungsi Tangga 7 knots 4.487083 3.431063Regresi Natural Cubic Spline
Secara Empiris
Menentukan basis Model Regresi Natural Cubic Spline.
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
df <- 6:12
set.seed(1501211024)
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
### menghitung rata-rata untuk 8 folds
mean_metric_spline3 <- colMeans(metric_spline3)
mean_metric_spline3
}
)
best_spline <- cbind(df=df,best_spline3)
# menampilkan hasil semua basis
best_spline## df rmse mae
## 1 6 4.361448 3.275680
## 2 7 4.332629 3.252226
## 3 8 4.348291 3.268900
## 4 9 4.358615 3.271120
## 5 10 4.329596 3.245079
## 6 11 4.366193 3.286472
## 7 12 4.368266 3.263797# menampilkan basis berdasarkan rsme dan mae yang terpilih
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)## df rmse mae
## 1 10 4.329596 3.245079
## 2 10 4.329596 3.245079Berdasarkan rsme dan mae, basis yang terpilih adalah df=10. Selanjutnya, menampilkan vektor knots dengan df=10.
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.7# menampilkan model
nilai.cv.ns <- rbind("NCS6" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("9 - NCS (df=10)")
model_ncs_h <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_h## Model rmse mae
## NCS6 9 - NCS (df=10) 4.329596 3.245079Secara Visualisasi
regnspline <- lm(mpg ~ ns(horsepower, knots=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7)))
prediksi <- predict(regnspline, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "NCS (df=10) horsepower vs mpg")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7),col="grey50",lty=2)
Perbandingan Model
Secara Empiris
rbind(model_poly_h, model_tangga_h, model_ncs_h)## Model rmse mae
## Poly2 Regresi Polinomial Derajat 2 4.382221 3.287275
## tangga Fungsi Tangga 7 knots 4.487083 3.431063
## NCS6 9 - NCS (df=10) 4.329596 3.245079Berdasarkan ketiga Model diatas, Model Natural Cubic Spline memiliki nilai rsme dan mae terkecil. Artinya Model Natural Cubic Spline adalah Model terbaik untuk data Auto dengan peubah horsepower vs mpg.
Secara Visualisasi
par(mfrow=c(2,2))
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Scatter Plot horsepower vs mpg")
polinom_d2 <- lm(mpg ~ poly(horsepower,2))
prediksi <- predict(polinom_d2, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 2")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
tangga8 <- lm(mpg ~ cut(horsepower,7))
prediksi <- predict(tangga8, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 11 knots")
ix <- sort(horsepower, index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l",lwd=2)
regnspline <- lm(mpg ~ ns(horsepower, knots=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7)))
prediksi <- predict(regnspline, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "NCS (df=10) horsepower vs mpg")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7),col="grey50",lty=2)
Pemodelan weight dan mpg
Menampilkan Data yang Digunakan
attach(Auto)
datatable(cbind(mpg, weight))plot(weight, mpg, pch=20, col="coral",
xlab = "Weight", ylab = "mile per gallon",
main = "Scatter Plot weight vs mpg")
Grafik di atas merupakan scatter plot dari 392 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah berat suatu kendaraan (Weight), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).
Secara grafik, terlihat bahwa pola data menunjukkan ke pola nonlinear. Untuk menemukan solusi optimal dari data kita bisa menggunakan Model Regresi Polinomial, Fungsi Tangga (Piecewise Constant) dan Natural Cubic Splines.
Regresi Polinomial
Secara Empiris
Cross Validation Data Auto menggunakan Model Regresi Polinomial
- Regresi Polinomial Derajat 2
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(weight,2,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly2)- Regresi Polinomial Derajat 3
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly3 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(weight,3,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly3)- Regresi Polinomial Derajat 4
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly4 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(weight,4,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly4)Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.
poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)
cvpoly <- rbind (poly2,poly3,poly4)
cvpoly## rmse mae
## poly2 4.183306 3.088252
## poly3 4.189345 3.093543
## poly4 4.201559 3.101442Berdasarkan hasil rata-rata rsme dan mae regresi polinomial diatas, terlihat bahwa Regresi Polinomial Derajat 2 memiliki nilai rsme dan mae paling kecil. Artinya, Regresi Polinomial Derajat 2 lebih baik dibandingkan Regresi Polinomial Derajat 3 dan 4.
# menampilkan model
nilai.cvpoly <- rbind("Poly2" = poly2)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 2")
model_poly_w <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_w## Model rmse mae
## Poly2 Regresi Polinomial Derajat 2 4.183306 3.088252Secara Visualisasi
par(mfrow=c(2,2))
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Scatter Plot weight vs mpg")
polinom_d2 <- lm(mpg ~ poly(weight,2))
prediksi <- predict(polinom_d2, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 2")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d3 <- lm(mpg ~ poly(weight,3))
prediksi <- predict(polinom_d3, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 3")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d4 <- lm(mpg ~ poly(weight,4))
prediksi <- predict(polinom_d4, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 4")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
Fungsi Tangga (Piecewise Constant)
Secara Empiris
Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
breaks <- 3:12
set.seed(1501211024)
best_tangga <- map_dfr(breaks, function(i){
metric_tangga <- 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_tangga
# menghitung rata-rata untuk 8 folds
mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)
# menampilkan semua jumlah interval
best_tangga## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.84 3.64
## 2 4.74 3.64
## 3 4.46 3.27
## 4 4.27 3.16
## 5 4.29 3.21
## 6 4.45 3.31
## 7 4.33 3.17
## 8 4.35 3.19
## 9 4.31 3.21
## 10 4.25 3.13# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)## breaks rmse mae
## 1 12 4.252092 3.130243
## 2 12 4.252092 3.130243Berdasarkan rmse dan mae jumlah interval terbaik adalah 12, maka pada kasus ini breaks 12 atau Fungsi Tangga dengan 11 knots adalah interval terbaik.
Secara Visualiasi
tangga12 <- lm(mpg ~ cut(weight,12))
prediksi <- predict(tangga12, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 11 knots")
ix <- sort(weight, index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
type="l",lwd=2)
# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 11 knots")
model_tangga_w <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_w ## Model rmse mae
## tangga Fungsi Tangga 11 knots 4.252092 3.130243Regresi Natural Cubic Spline
Secara Empiris
Menentukan basis Model Regresi Natural Cubic Spline.
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
df <- 6:12
set.seed(1501211024)
best_spline3 <- map_dfr(df, function(i){
metric_spline3 <- 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_spline3
### menghitung rata-rata untuk 8 folds
mean_metric_spline3 <- colMeans(metric_spline3)
mean_metric_spline3
}
)
best_spline <- cbind(df=df,best_spline3)
# menampilkan hasil semua basis
best_spline## df rmse mae
## 1 6 4.198793 3.119335
## 2 7 4.167602 3.072446
## 3 8 4.184381 3.095133
## 4 9 4.206877 3.096451
## 5 10 4.187345 3.084281
## 6 11 4.176897 3.085570
## 7 12 4.180123 3.096528# menampilkan basis berdasarkan rsme dan mae yang terpilih
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)## df rmse mae
## 1 7 4.167602 3.072446
## 2 7 4.167602 3.072446Berdasarkan rsme dan mae, basis yang terpilih adalah df=7. Selanjutnya, menampilkan vektor knots dengan df=7.
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.286# menampilkan model
nilai.cv.ns <- rbind("NCS6" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("6 - NCS (df=7)")
model_ncs_w <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_w## Model rmse mae
## NCS6 6 - NCS (df=7) 4.167602 3.072446Secara Visualisasi
regnspline <- lm(mpg ~ ns(weight, knots=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "NCS (df=7) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286),col="grey50",lty=2)
Perbandingan Model
Secara Empiris
rbind(model_poly_w, model_tangga_w, model_ncs_w)## Model rmse mae
## Poly2 Regresi Polinomial Derajat 2 4.183306 3.088252
## tangga Fungsi Tangga 11 knots 4.252092 3.130243
## NCS6 6 - NCS (df=7) 4.167602 3.072446Berdasarkan ketiga Model diatas, Model Natural Cubic Spline memiliki nilai rsme dan mae terkecil. Artinya Model Natural Cubic Spline adalah Model terbaik untuk data Auto dengan peubah weight vs mpg.
Secara Visualisasi
par(mfrow=c(2,2))
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Scatter Plot weight vs mpg")
polinom_d2 <- lm(mpg ~ poly(weight,2))
prediksi <- predict(polinom_d2, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 2")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
tangga12 <- lm(mpg ~ cut(weight,12))
prediksi <- predict(tangga12, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 11 knots")
ix <- sort(weight, index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
type="l",lwd=2)
regnspline <- lm(mpg ~ ns(weight, knots=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "NCS (df=7) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286),col="grey50",lty=2)
Pemodelan displacement dan mpg
Menampilkan Data yang Digunakan
attach(Auto)
datatable(cbind(mpg, displacement))plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Scatter Plot displacement vs mpg")
Grafik di atas merupakan scatter plot dari 392 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah jarak tempuh suatu kendaraan (displacement), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).
Secara grafik, terlihat bahwa pola data menunjukkan ke pola nonlinear. Untuk menemukan solusi optimal dari data kita bisa menggunakan Model Regresi Polinomial, Fungsi Tangga (Piecewise Constant) dan Natural Cubic Splines.
Regresi Polinomial
Secara Empiris
Cross Validation Data Auto menggunakan Model Regresi Polinomial
- Regresi Polinomial Derajat 2
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(displacement,2,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly2)- Regresi Polinomial Derajat 3
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly3 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(displacement,3,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly3)- Regresi Polinomial Derajat 4
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
set.seed(1501211024)
metric_poly4 <- map_dfr(cross_val$splits,
function(x){
#data training
mod <- lm(mpg ~ poly(displacement,4,raw = T),
data=Auto[x$in_id,])
#data testing
pred <- predict(mod, newdata=Auto[-x$in_id,])
#nilai y pada data testing
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)
}
)
datatable(metric_poly4)Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.
poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)
cvpoly <- rbind (poly2,poly3,poly4)
cvpoly## rmse mae
## poly2 4.331144 3.176381
## poly3 4.330340 3.176739
## poly4 4.332225 3.176730Berdasarkan hasil rata-rata rsme dan mae regresi polinomial diatas, terlihat bahwa Regresi Polinomial Derajat 3 memiliki nilai rsme dan mae paling kecil. Artinya, Regresi Polinomial Derajat 3 lebih baik dibandingkan Regresi Polinomial Derajat 2 dan 4.
# menampilkan model
nilai.cvpoly <- rbind("Poly3" = poly3)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 3")
model_poly_d <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_d## Model rmse mae
## Poly3 Regresi Polinomial Derajat 3 4.33034 3.176739Secara Visualisasi
par(mfrow=c(2,2))
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Scatter Plot displacement vs mpg")
polinom_d2 <- lm(mpg ~ poly(displacement,2))
prediksi <- predict(polinom_d2, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 2")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d3 <- lm(mpg ~ poly(displacement,3))
prediksi <- predict(polinom_d3, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 3")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
polinom_d4 <- lm(mpg ~ poly(displacement,4))
prediksi <- predict(polinom_d4, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 4")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
Fungsi Tangga (Piecewise Constant)
Secara Empiris
Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
breaks <- 3:12
set.seed(1501211024)
best_tangga <- map_dfr(breaks, function(i){
metric_tangga <- 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_tangga
# menghitung rata-rata untuk 8 folds
mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)
# menampilkan semua jumlah interval
best_tangga## # A tibble: 10 x 2
## rmse mae
## <dbl> <dbl>
## 1 4.91 3.69
## 2 4.81 3.64
## 3 4.69 3.52
## 4 4.60 3.41
## 5 4.58 3.39
## 6 4.39 3.25
## 7 4.23 3.10
## 8 4.24 3.10
## 9 4.37 3.25
## 10 4.42 3.27# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)## breaks rmse mae
## 1 9 4.234784 3.099262
## 2 9 4.234784 3.099262Berdasarkan rmse dan mae jumlah interval terbaik adalah 9, maka pada kasus ini breaks 9 atau Fungsi Tangga dengan 8 knots adalah interval terbaik.
Secara Visualiasi
tangga9 <- lm(mpg ~ cut(displacement,8))
prediksi <- predict(tangga9, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 8 knots")
ix <- sort(displacement, index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l",lwd=2)
# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 8 knots")
model_tangga_d <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_d ## Model rmse mae
## tangga Fungsi Tangga 8 knots 4.234784 3.099262Regresi Natural Cubic Spline
Secara Empiris
Menentukan basis Model Regresi Natural Cubic Spline.
set.seed(1501211024)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")
df <- 6:12
set.seed(1501211024)
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
### menghitung rata-rata untuk 8 folds
mean_metric_spline3 <- colMeans(metric_spline3)
mean_metric_spline3
}
)
best_spline <- cbind(df=df,best_spline3)
# menampilkan hasil semua basis
best_spline## df rmse mae
## 1 6 4.231043 3.137840
## 2 7 4.187671 3.119895
## 3 8 4.192260 3.147345
## 4 9 4.195433 3.155416
## 5 10 4.210086 3.170656
## 6 11 4.174756 3.109714
## 7 12 4.176774 3.089360# menampilkan basis berdasarkan rsme dan mae yang terpilih
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)## df rmse mae
## 1 11 4.174756 3.109714
## 2 12 4.176774 3.089360Berdasarkan rsme basis yang terpilih adalah df=11, dan berdasarkan mae basis yang terpilih adalah df=12. Namun pada kasus ini basis yang dipilih berdasarkan rsme terkecil yaitu df=11. Selanjutnya, menampilkan vektor knots dengan df=11 dan 12.
attr(ns(Auto$displacement, df=11), "knots")## 9.090909% 18.18182% 27.27273% 36.36364% 45.45455% 54.54545% 63.63636% 72.72727%
## 90 97 107 120 140 168 231 258
## 81.81818% 90.90909%
## 318 351# menampilkan model
nilai.cv.ns <- rbind("NCS10" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("10 - NCS (df=11)")
model_ncs_d <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_d## Model rmse mae
## NCS10 10 - NCS (df=11) 4.174756 3.109714Secara Visualisasi
regnspline <- lm(mpg ~ ns(displacement, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "NCS (df=11) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)
Perbandingan Model
Secara Empiris
rbind(model_poly_d, model_tangga_d, model_ncs_d)## Model rmse mae
## Poly3 Regresi Polinomial Derajat 3 4.330340 3.176739
## tangga Fungsi Tangga 8 knots 4.234784 3.099262
## NCS10 10 - NCS (df=11) 4.174756 3.109714Secara empiris jika melihat nilai rsme dan mae, Model Natural Cubic Spline memiliki nilai rsme terkecil, Fungsi Tangga 8 knots memiliki nilai mae terkecil. Jika kita memilih Model terbaik berdasarkan nilai rsme terkecil maka Model Natural Cubic Spline adalah Model terbaik untuk data Auto untuk peubah displacement vs mpg.
Secara Visualisasi
par(mfrow=c(2,2))
plot(weight, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Scatter Plot displacement vs mpg")
polinom_d3 <- lm(mpg ~ poly(displacement,3))
prediksi <- predict(polinom_d3, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Regresi Polinomial Ordo 3")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
main = "Polynomial Regression", col="blue", type="l",lwd=2)
tangga9 <- lm(mpg ~ cut(displacement,8))
prediksi <- predict(tangga9, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 8 knots")
ix <- sort(displacement, index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l",lwd=2)
regnspline <- lm(mpg ~ ns(displacement, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "NCS (df=11) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)
Sedangkan secara visualisasi, Model Fungsi Tangga 8 knots memiliki grafik yang lebih baik dibandingkan Model lainnya.
Perbandingan semua model dan semua data
Secara Empiris
horsepower1 <- rbind(horsepower1 = rbind(model_poly_h, model_tangga_h, model_ncs_h) %>% slice_min(rmse) )
weight1 <- rbind(weght1 = rbind(model_poly_w, model_tangga_w, model_ncs_w) %>% slice_min(rmse))
displacement1 <- rbind(displacement1 = rbind(model_poly_d, model_tangga_d, model_ncs_d) %>% slice_min(rmse))
rbind(horsepower1,weight1,displacement1)## Model rmse mae
## horsepower1 9 - NCS (df=10) 4.329596 3.245079
## weght1 6 - NCS (df=7) 4.167602 3.072446
## displacement1 10 - NCS (df=11) 4.174756 3.109714Bedasarkan perbandingan model dengan melihat rsme, Model Natural Cubic Spline adalah model terbaik untuk pada data Auto dengan peubah horsepower vs mpg, weight vs mpg, dan displacement vs mpg.
Secara Visualisasi
par(mfrow=c(2,2))
regnspline <- lm(mpg ~ ns(horsepower, knots=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7)))
prediksi <- predict(regnspline, data.frame(horsepower), interval = "predict")
plot(horsepower, mpg, pch=20, col="coral",
xlab = "horsepower", ylab = "mile per gallon",
main = "NCS (df=10) horsepower vs mpg")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7),col="grey50",lty=2)
regnspline <- lm(mpg ~ ns(weight, knots=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
xlab = "weight", ylab = "mile per gallon",
main = "NCS (df=7) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286),col="grey50",lty=2)
tangga9 <- lm(mpg ~ cut(displacement,8))
prediksi <- predict(tangga9, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "Fungsi Tangga dengan 8 knots")
ix <- sort(displacement, index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l",lwd=2)
regnspline <- lm(mpg ~ ns(displacement, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
xlab = "displacement", ylab = "mile per gallon",
main = "NCS (df=11) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)