Data disimulasikan dengan hubungan non-linier (kuadratik) antara
x dan y.
set.seed(123)
n <- 100
x <- seq(0, 10, length.out = n)
y <- 5 + 2*x - 0.3*x^2 + rnorm(n, mean = 0, sd = 2)
data <- data.frame(x = x, y = y)
head(data)
## x y
## 1 0.0000000 3.879049
## 2 0.1010101 4.738604
## 3 0.2020202 8.509213
## 4 0.3030303 5.719529
## 5 0.4040404 6.017682
## 6 0.5050505 9.363708
#Interpretasi
Berdasarkan output diatas, diperoleh 100 data bangkitan untuk variabel x dan y, dengan variabel x memiliki jarak yang sama dimulai dari 0 hingga 10. Kemudian, didapatkan variabel y berdistribusi normal dengan rata-rata 0 dan standar deviasi 2.
Tiga model dicocokkan sebagai pembanding: linier, polinomial derajat 2, dan derajat 3.
model_linier <- lm(y ~ x, data = data)
summary(model_linier)
##
## Call:
## lm(formula = y ~ x, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.4317 -1.4051 -0.0466 2.2166 6.6109
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.88171 0.55369 17.847 <2e-16 ***
## x -0.95028 0.09566 -9.934 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.789 on 98 degrees of freedom
## Multiple R-squared: 0.5017, Adjusted R-squared: 0.4967
## F-statistic: 98.68 on 1 and 98 DF, p-value: < 2.2e-16
Berdasarkan model tersebut, konstanta sebesar 9,88171 menyatakan bahwa tanpa dipengaruhi oleh variabel x maka variabel y adalah sebesar 9,88171 satuan. Koefisien regresi sebesar -0,95028 menyatakan bahwa setiap penambahan 1 satuan variabel x akan menurunkan variabel y sebesar 0,95028.
Koefisien determinasi sebesar 0,5017 menyatakan bahwa sebesar 50,17% variasi yang terjadi pada variabel y disebabkan oleh variabel x. Hal ini menandakan bahwa ada sekitar 49,83% variasi yang terjadi disebabkan oleh variabel lain yang tidak diketahui.
model_poli2 <- lm(y ~ poly(x, 2, raw = TRUE), data = data)
summary(model_poli2)
##
## Call:
## lm(formula = y ~ poly(x, 2, raw = TRUE), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8136 -1.1977 -0.0533 1.3549 4.3891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.33982 0.53777 9.929 <2e-16 ***
## poly(x, 2, raw = TRUE)1 1.80266 0.24855 7.253 1e-10 ***
## poly(x, 2, raw = TRUE)2 -0.27529 0.02405 -11.447 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.829 on 97 degrees of freedom
## Multiple R-squared: 0.788, Adjusted R-squared: 0.7837
## F-statistic: 180.3 on 2 and 97 DF, p-value: < 2.2e-16
Berdasarkan model polinomial 2 diatas, konstanta sebesar 5,33982 menyatakan bahwa tanpa dipengaruhi oleh variabel x maka variabel y adalah sebesar 5,33982 satuan. Koefisien regresi sebesar 1,80266 menyatakan bahwa setiap penambahan 1 satuan variabel x akan menaikkan variabel y sebesar 1,80266. Koefisien regresi sebesar -0,27529 menyatakan bahwa setiap penambahan 1 satuan variabel x² akan menurunkan variabel y sebesar 0,27529
Koefisien determinasi sebesar 0,788 menyatakan bahwa sebesar 78,8% variasi yang terjadi pada variabel y disebabkan oleh variabel x dan x². Hal ini menandakan bahwa ada sekitar 21,2% variasi yang terjadi disebabkan oleh variabel lain yang tidak diketahui.
model_poli3 <- lm(y ~ poly(x, 3, raw = TRUE), data = data)
summary(model_poli3)
##
## Call:
## lm(formula = y ~ poly(x, 3, raw = TRUE), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7866 -1.1996 -0.0497 1.3386 4.3777
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.282801 0.708433 7.457 3.93e-11 ***
## poly(x, 3, raw = TRUE)1 1.872854 0.616613 3.037 0.00307 **
## poly(x, 3, raw = TRUE)2 -0.292931 0.143693 -2.039 0.04424 *
## poly(x, 3, raw = TRUE)3 0.001176 0.009443 0.125 0.90117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.838 on 96 degrees of freedom
## Multiple R-squared: 0.7881, Adjusted R-squared: 0.7815
## F-statistic: 119 on 3 and 96 DF, p-value: < 2.2e-16
Berdasarkan model polinomial 3 diatas, konstanta sebesar 5,282801 menyatakan bahwa tanpa dipengaruhi oleh variabel x maka variabel y adalah sebesar 5,282801 satuan. Koefisien regresi sebesar 1,872854 menyatakan bahwa setiap penambahan 1 satuan variabel x akan menaikkan variabel y sebesar 1,872854. Koefisien regresi sebesar -0,292931 menyatakan bahwa setiap penambahan 1 satuan variabel x² akan menurunkan variabel y sebesar 0,292931. Koefisien regresi sebesar 0,001176 menyatakan bahwa setiap penambahan 1 satuan variabel x³ akan menaikkan variabel y sebesar 0,001176.
Koefisien determinasi sebesar 0,7881 menyatakan bahwa sebesar 78,81% variasi yang terjadi pada variabel y disebabkan oleh variabel x, x² dan x³. Hal ini menandakan bahwa ada sekitar 21,19% variasi yang terjadi disebabkan oleh variabel lain yang tidak diketahui.
anova(model_linier, model_poli2, model_poli3)
## Analysis of Variance Table
##
## Model 1: y ~ x
## Model 2: y ~ poly(x, 2, raw = TRUE)
## Model 3: y ~ poly(x, 3, raw = TRUE)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 98 762.42
## 2 97 324.33 1 438.09 129.6931 <2e-16 ***
## 3 96 324.28 1 0.05 0.0155 0.9012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Interpretasi
Berdasarkan hasil uji ANOVA yang membandingkan model linier, polinomial derajat dua, dan polinomial derajat tiga, diketahui bahwa model linier menghasilkan RSS sebesar 762,42. Penambahan suku kuadratik pada model polinomial derajat dua menurunkan RSS secara menjadi 324,33, dengan selisih penurunan sebesar 438,09 dan nilai p-value < 2e-16 yang mengindikasikan signifikansi statistik yang sangat kuat. Hal ini membuktikan bahwa model polinomial derajat dua jauh lebih baik dalam menjelaskan variabilitas data dibandingkan model linier. Sebaliknya, peningkatan ke polinomial derajat tiga hanya menurunkan RSS menjadi 324,28, dengan penurunan marginal sebesar 0,05 dan nilai p = 0,9012 yang tidak signifikan. Dengan demikian, model polinomial derajat dua dipilih sebagai model terbaik karena memberikan peningkatan performa yang signifikan tanpa kompleksitas tambahan yang tidak perlu. Model akhir yang diperoleh adalah ŷ = 5,33982 + 1,80266x − 0,27529x².
AIC(model_linier, model_poli2, model_poli3)
## df AIC
## model_linier 3 492.9204
## model_poli2 4 409.4469
## model_poli3 5 411.4307
#Interpretasi
Berdasarkan output diatas, dapat diketahui bahwa model linear memiliki nilai AIC paling besar dengan nilai sebesar 492,9204. Kemudian dilanjutkan dengan model polinomial derajat 3 dengan nilai AIC sebesar 411,4307, dan model yang memiliki nilai AIC paling kecil adalah model polinomial derajat 2 dengan nilai AIC sebesar 409,4469.
data$pred_linier <- predict(model_linier)
data$pred_poli2 <- predict(model_poli2)
data$pred_poli3 <- predict(model_poli3)
ggplot(data, aes(x = x, y = y)) +
geom_point(alpha = 0.5) +
geom_line(aes(y = pred_linier, color = "Linier"), linewidth = 1) +
geom_line(aes(y = pred_poli2, color = "Polinomial derajat 2"), linewidth = 1) +
geom_line(aes(y = pred_poli3, color = "Polinomial derajat 3"), linewidth = 1) +
labs(title = "Perbandingan Regresi Linier vs Polinomial",
x = "X", y = "Y", color = "Model") +
theme_minimal()
#Interpretasi
Berdasarkan visualisasi di atas, terlihat bahwa model linier (ditandai dengan garis merah) membentuk pola garis lurus dengan hubungan negatif atau menurun. Sementara itu, model polinomial derajat 2 (biru) dan derajat 3 (hijau) keduanya membentuk kurva melengkung. Kedua kurva tersebut hampir saling menutupi karena pola yang dihasilkan nyaris identik, menunjukkan bahwa penambahan suku pangkat tiga tidak memberikan perubahan yang berarti pada bentuk kurva dibandingkan model derajat dua.
rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2))
data.frame(
Model = c("Linier", "Polinomial derajat 2", "Polinomial derajat 3"),
RMSE = c(rmse(data$y, data$pred_linier),
rmse(data$y, data$pred_poli2),
rmse(data$y, data$pred_poli3))
)
## Model RMSE
## 1 Linier 2.761195
## 2 Polinomial derajat 2 1.800917
## 3 Polinomial derajat 3 1.800772
#Interpretasi
Berdasarkan hasil perhitungan RMSE, model linier memiliki nilai kesalahan tertinggi yaitu sebesar 2,761195, yang menunjukkan performa prediksi paling rendah dibandingkan kedua model lainnya. Sementara itu, model polinomial derajat 2 dan derajat 3 memiliki nilai RMSE yang hampir sama, masing-masing sebesar 1,800917 dan 1,800772. Kedekatan nilai ini mengindikasikan bahwa kedua model memiliki tingkat akurasi yang setara. Dengan demikian, hasil ini memperkuat kesimpulan sebelumnya bahwa model polinomial derajat 2 merupakan model yang optimal, karena memberikan performa yang sama baiknya dengan model derajat 3 namun dengan kompleksitas yang lebih rendah.
Skema 10-fold cross-validation diulang 5 kali untuk menguji derajat polinomial 1 sampai 6.
names(data)
## [1] "x" "y" "pred_linier" "pred_poli2" "pred_poli3"
library(caret)
set.seed(42)
kontrol <- trainControl(
method = "repeatedcv",
number = 10,
repeats = 5
)
derajat_max <- 6
hasil_cv <- data.frame(
derajat = integer(),
RMSE = numeric()
)
for (d in 1:derajat_max) {
data_cv <- data
# Membuat variabel polynomial
if (d == 1) {
data_cv$x_poly <- data_cv$x
} else {
poly_x <- poly(data_cv$x, degree = d, raw = TRUE)
data_cv$x_poly <- poly_x[, 1]
for (j in 2:d) {
data_cv[[paste0("x_poly", j)]] <- poly_x[, j]
}
}
# Formula sesuai derajat
if (d == 1) {
formula_model <- y ~ x_poly
} else {
variabel <- paste0("x_poly", 2:d)
formula_model <- as.formula(
paste("y ~ x_poly +", paste(variabel, collapse = " + "))
)
}
model_cv <- train(
formula_model,
data = data_cv,
method = "lm",
trControl = kontrol
)
hasil_cv <- rbind(
hasil_cv,
data.frame(
derajat = d,
RMSE = min(model_cv$results$RMSE)
)
)
}
hasil_cv
## derajat RMSE
## 1 1 2.757805
## 2 2 1.825068
## 3 3 1.859968
## 4 4 1.877628
## 5 5 1.871701
## 6 6 1.880719
derajat_optimal <- hasil_cv$derajat[which.min(hasil_cv$RMSE)]
cat("Derajat optimal berdasarkan CV:", derajat_optimal, "\n")
## Derajat optimal berdasarkan CV: 2
ggplot(hasil_cv, aes(x = derajat, y = RMSE)) +
geom_line(
color = "steelblue",
linewidth = 1
) +
geom_point(
size = 3,
color = "steelblue"
) +
geom_vline(
xintercept = derajat_optimal,
linetype = "dashed",
color = "firebrick"
) +
scale_x_continuous(
breaks = 1:derajat_max
) +
labs(
title = "Pemilihan Derajat Optimal via 10-Fold Cross-Validation",
x = "Derajat Polinomial",
y = "RMSE Rata-rata (Validasi Silang)"
) +
theme_minimal()
#Interpretasi
Berdasarkan grafik, RMSE menurun dari derajat 1 ke derajat 2, kemudian meningkat lagi pada derajat 3 hingga 6. Hal ini menunjukkan bahwa derajat 3 memiliki RMSE lebih besar daripada derajat 2, sehingga derajat 2 merupakan model yang optimal.
model_final <- lm(y ~ poly(x, derajat_optimal, raw = TRUE), data = data)
summary(model_final)
##
## Call:
## lm(formula = y ~ poly(x, derajat_optimal, raw = TRUE), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8136 -1.1977 -0.0533 1.3549 4.3891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.33982 0.53777 9.929 <2e-16 ***
## poly(x, derajat_optimal, raw = TRUE)1 1.80266 0.24855 7.253 1e-10 ***
## poly(x, derajat_optimal, raw = TRUE)2 -0.27529 0.02405 -11.447 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.829 on 97 degrees of freedom
## Multiple R-squared: 0.788, Adjusted R-squared: 0.7837
## F-statistic: 180.3 on 2 and 97 DF, p-value: < 2.2e-16
Berdasarkan output diatas, dapat diketahui bahwa model polinomial derajat 2 terpilih sebagai model terbaik. Hal ini didukung oleh nilai R-squared sebesar 0,788 yang mengindikasikan kemampuan prediksi yang baik, koefisien yang signifikan secara statistik, serta nilai AIC dan RMSE yang paling rendah dibandingkan model lainnya. Kombinasi ini menjadikan model derajat 2 sebagai pilihan yang paling optimal.
#Memunculkan nilai Y topi (Y prediksi)
model_final <- lm(y ~ poly(x, 2, raw = TRUE), data = data)
#Hitung nilai prediksi y
data$Y_topi <- predict(model_final)
#tampilkan kolom x,y dan y_topi
i <- data[, c("x", "y", "Y_topi")] # <-- MASALAH ADA DI SINI
head(i)
## x y Y_topi
## 1 0.0000000 3.879049 5.339823
## 2 0.1010101 4.738604 5.519101
## 3 0.2020202 8.509213 5.692761
## 4 0.3030303 5.719529 5.860804
## 5 0.4040404 6.017682 6.023229
## 6 0.5050505 9.363708 6.180036
Residual <- residuals(model_final)
# Atau tampilkan dalam dataframe
data$Residual <- residuals(model_final)
print(data[, c("x", "y", "Y_topi", "Residual")])
## x y Y_topi Residual
## 1 0.0000000 3.87904871 5.33982287 -1.460774162
## 2 0.1010101 4.73860431 5.51910079 -0.780496476
## 3 0.2020202 8.50921338 5.69276104 2.816452348
## 4 0.3030303 5.71952918 5.86080361 -0.141274434
## 5 0.4040404 6.01768168 6.02322852 -0.005546837
## 6 0.5050505 9.36370818 6.18003576 3.183672422
## 7 0.6060606 7.02376079 6.33122532 0.692535462
## 8 0.7070707 3.73403425 6.47679722 -2.742762972
## 9 0.8080808 5.04655753 6.61675145 -1.570193912
## 10 0.9090909 5.67892399 6.75108800 -1.072164008
## 11 1.0101010 9.16227440 6.87980689 2.282467513
## 12 1.1111111 7.57147951 7.00290810 0.568571404
## 13 1.2121212 7.78501398 7.12039165 0.664622330
## 14 1.3131313 7.33033390 7.23225752 0.098076385
## 15 1.4141414 6.11666178 7.33850572 -1.221843946
## 16 1.5151515 10.91542407 7.43913626 3.476287814
## 17 1.6161616 8.44443068 7.53414912 0.910281560
## 18 1.7171717 3.61650551 7.62354431 -4.007038800
## 19 1.8181818 9.04733990 7.70732183 1.340018071
## 20 1.9191919 6.78781174 7.78548168 -0.997669946
## 21 2.0202020 5.68039177 7.85802386 -2.177632095
## 22 2.1212121 7.45661215 7.92494837 -0.468336219
## 23 2.2222222 5.91095407 7.98625521 -2.075301147
## 24 2.3232323 6.56945966 8.04194438 -1.472484723
## 25 2.4242424 6.83532091 8.09201588 -1.256694968
## 26 2.5252525 4.76404833 8.13646971 -3.372421376
## 27 2.6262626 9.85892273 8.17530587 1.683616859
## 28 2.7272727 8.52988673 8.20852436 0.321362376
## 29 2.8282828 5.98053666 8.23612517 -2.255588517
## 30 2.9292929 10.79198858 8.25810832 2.533880261
## 31 3.0303030 9.15871357 8.27447380 0.884239771
## 32 3.1313131 7.73094672 8.28522160 -0.554274883
## 33 3.2323232 10.12052374 8.29035174 1.830172006
## 34 3.3333333 10.08960031 8.28986420 1.799736106
## 35 3.4343434 9.97343458 8.28375900 1.689675587
## 36 3.5353535 9.69837019 8.27203612 1.426334071
## 37 3.6363636 9.41362043 8.25469558 1.158924856
## 38 3.7373737 8.16053532 8.23173736 -0.071202041
## 39 3.8383838 7.64488520 8.20316147 -0.558276270
## 40 3.9393939 7.46219849 8.16896791 -0.706769420
## 41 4.0404040 6.79393468 8.12915669 -1.335222006
## 42 4.1414141 7.72160040 8.08372779 -0.362127388
## 43 4.2424242 5.55460675 8.03268122 -2.478074473
## 44 4.3434343 12.36515405 7.97601698 4.389137069
## 45 4.4444444 10.37888696 7.91373507 2.465151889
## 46 4.5454545 5.64634482 7.84583549 -2.199490673
## 47 4.6464646 7.01026951 7.77231824 -0.762048731
## 48 4.7474747 6.80008384 7.69318332 -0.893099475
## 49 4.8484848 9.20455834 7.60843073 1.596127608
## 50 4.9494949 7.38300169 7.51806047 -0.135058776
## 51 5.0505051 7.95536675 7.42207253 0.533294216
## 52 5.1515152 7.28450429 7.32046693 -0.035962646
## 53 5.2525253 7.14260313 7.21324366 -0.070640526
## 54 5.3535354 9.84617304 7.10040271 2.745770326
## 55 5.4545455 6.53192910 6.98194410 -0.450014998
## 56 5.5555556 9.88479306 6.85786782 3.026925244
## 57 5.6565657 3.61660520 6.72817386 -3.111568665
## 58 5.7575758 7.73947543 6.59286224 1.146613197
## 59 5.8585859 6.66797173 6.45193294 0.216038785
## 60 5.9595960 6.69603986 6.30538598 0.390653881
## 61 6.0606061 6.86120734 6.15322134 0.707986002
## 62 6.1616162 4.92893130 5.99543903 -1.066507732
## 63 6.2626263 5.09269145 5.83203905 -0.739347608
## 64 6.3636364 3.54136163 5.66302141 -2.121659776
## 65 6.4646465 3.24821430 5.48838609 -2.240171786
## 66 6.5656566 5.80601657 5.30813310 0.497883473
## 67 6.6666667 5.89641956 5.12226244 0.774157116
## 68 6.7676768 4.90092734 4.93077411 -0.029846772
## 69 6.8686869 6.42825088 4.73366811 1.694582772
## 70 6.9696970 8.46656056 4.53094444 3.935616115
## 71 7.0707071 3.16088227 4.32260310 -1.161720835
## 72 7.1717172 -0.70496157 4.10864409 -4.813605655
## 73 7.2727273 5.68916300 3.88906741 1.800095591
## 74 7.3737374 2.01747237 3.66387306 -1.646400691
## 75 7.4747475 1.81192277 3.43306103 -1.621138260
## 76 7.5757576 4.98502704 3.19663134 1.788395696
## 77 7.6767677 2.10416075 2.95458398 -0.850423227
## 78 7.7777778 -0.03402802 2.70691894 -2.740946961
## 79 7.8787879 2.49759319 2.45363624 0.043956946
## 80 7.9797980 1.57866050 2.19473586 -0.616075370
## 81 8.0808081 1.58330676 1.93021782 -0.346911058
## 82 8.1818182 2.05155254 1.66008210 0.391470433
## 83 8.2828283 0.24276319 1.38432872 -1.141565525
## 84 8.3838384 1.96980605 1.10295766 0.866848389
## 85 8.4848485 -0.06907230 0.81596894 -0.885041233
## 86 8.5858586 0.72019080 0.52336254 0.196828264
## 87 8.6868687 1.92890913 0.22513847 1.703770654
## 88 8.7878788 0.27807648 -0.07870327 0.356779747
## 89 8.8888889 -1.57778910 -0.38816268 -1.189626421
## 90 8.9898990 1.03192806 -0.70323975 1.735167817
## 91 9.0909091 0.37543746 -1.02393450 1.399371968
## 92 9.1919192 -0.86678123 -1.35024692 0.483465697
## 93 9.2929293 -1.84423840 -1.68217701 -0.162061383
## 94 9.3939394 -3.94176257 -2.01972477 -1.922037791
## 95 9.4949495 -0.33501589 -2.36289021 2.027874319
## 96 9.5959596 -4.63333215 -2.71167331 -1.921658846
## 97 9.6969697 0.55923899 -3.06607408 3.625313067
## 98 9.7979798 -1.13894159 -3.42609252 2.287150932
## 99 9.8989899 -5.07042123 -3.79172863 -1.278692593
## 100 10.0000000 -7.05284180 -4.16298242 -2.889859384
Berdasarkan hasil analisis regresi polinomial pada 100 data observasi, nilai prediksi (ŷ) menunjukkan pola hubungan berbentuk parabola terbalik antara variabel x dan y. Nilai y diperkirakan meningkat seiring bertambahnya x hingga mencapai titik puncak sekitar x = 3,33 dengan prediksi ŷ = 8,29, kemudian menurun secara bertahap hingga titik terendah di x = 10 dengan prediksi ŷ = -4,16. Pola ini mengonfirmasi bahwa model mampu menangkap hubungan non-linear yang sesuai dengan persamaan teoritis y = 5 + 2x - 0,3x²
poly(x, derajat, raw = TRUE) digunakan agar koefisien
dapat diinterpretasi langsung sebagai \(\beta_1 x + \beta_2 x^2 + \dots\)x distandardisasi terlebih dahulu
(scale(x)) agar model lebih stabil secara numerik.