Regresi Polinomial dengan R
Contoh Penerapan dan Pemilihan Derajat via Cross-Validation
Rini Rahmah
2026-09-02
- Simulasi Data
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- Pencocokan Model
Tiga model dicocokkan sebagai pembanding: linier, polinomial derajat 2, dan derajat 3.
model_linier <- lm(y ~ x, data = data)
model_poli2 <- lm(y ~ poly(x, 2, raw = TRUE), data = data)
model_poli3 <- lm(y ~ poly(x, 3, raw = TRUE), data = data)
# y = nilai aktual (hasil simulasi)
# ŷ (y topi) = nilai prediksi/fitted value, akan ditambahkan setelah model dicocokkan
data$y_hat_linier <- predict(model_linier)
data$y_hat_linier## [1] 9.8817094 9.7857215 9.6897335 9.5937456 9.4977577 9.4017697 9.3057818
## [8] 9.2097939 9.1138060 9.0178180 8.9218301 8.8258422 8.7298542 8.6338663
## [15] 8.5378784 8.4418904 8.3459025 8.2499146 8.1539266 8.0579387 7.9619508
## [22] 7.8659628 7.7699749 7.6739870 7.5779990 7.4820111 7.3860232 7.2900352
## [29] 7.1940473 7.0980594 7.0020714 6.9060835 6.8100956 6.7141076 6.6181197
## [36] 6.5221318 6.4261439 6.3301559 6.2341680 6.1381801 6.0421921 5.9462042
## [43] 5.8502163 5.7542283 5.6582404 5.5622525 5.4662645 5.3702766 5.2742887
## [50] 5.1783007 5.0823128 4.9863249 4.8903369 4.7943490 4.6983611 4.6023731
## [57] 4.5063852 4.4103973 4.3144093 4.2184214 4.1224335 4.0264455 3.9304576
## [64] 3.8344697 3.7384818 3.6424938 3.5465059 3.4505180 3.3545300 3.2585421
## [71] 3.1625542 3.0665662 2.9705783 2.8745904 2.7786024 2.6826145 2.5866266
## [78] 2.4906386 2.3946507 2.2986628 2.2026748 2.1066869 2.0106990 1.9147110
## [85] 1.8187231 1.7227352 1.6267472 1.5307593 1.4347714 1.3387834 1.2427955
## [92] 1.1468076 1.0508196 0.9548317 0.8588438 0.7628559 0.6668679 0.5708800
## [99] 0.4748921 0.3789041data$y_hat_poli2 <- predict(model_poli2)
data$y_hat_poli2## [1] 5.33982287 5.51910079 5.69276104 5.86080361 6.02322852 6.18003576
## [7] 6.33122532 6.47679722 6.61675145 6.75108800 6.87980689 7.00290810
## [13] 7.12039165 7.23225752 7.33850572 7.43913626 7.53414912 7.62354431
## [19] 7.70732183 7.78548168 7.85802386 7.92494837 7.98625521 8.04194438
## [25] 8.09201588 8.13646971 8.17530587 8.20852436 8.23612517 8.25810832
## [31] 8.27447380 8.28522160 8.29035174 8.28986420 8.28375900 8.27203612
## [37] 8.25469558 8.23173736 8.20316147 8.16896791 8.12915669 8.08372779
## [43] 8.03268122 7.97601698 7.91373507 7.84583549 7.77231824 7.69318332
## [49] 7.60843073 7.51806047 7.42207253 7.32046693 7.21324366 7.10040271
## [55] 6.98194410 6.85786782 6.72817386 6.59286224 6.45193294 6.30538598
## [61] 6.15322134 5.99543903 5.83203905 5.66302141 5.48838609 5.30813310
## [67] 5.12226244 4.93077411 4.73366811 4.53094444 4.32260310 4.10864409
## [73] 3.88906741 3.66387306 3.43306103 3.19663134 2.95458398 2.70691894
## [79] 2.45363624 2.19473586 1.93021782 1.66008210 1.38432872 1.10295766
## [85] 0.81596894 0.52336254 0.22513847 -0.07870327 -0.38816268 -0.70323975
## [91] -1.02393450 -1.35024692 -1.68217701 -2.01972477 -2.36289021 -2.71167331
## [97] -3.06607408 -3.42609252 -3.79172863 -4.16298242data$y_hat_poli3 <- predict(model_poli3)
data$y_hat_poli3## [1] 5.28280126 5.46899089 5.64921021 5.82346649 5.99176701 6.15411903
## [7] 6.31052982 6.46100666 6.60555682 6.74418756 6.87690616 7.00371990
## [13] 7.12463603 7.23966183 7.34880458 7.45207154 7.54946998 7.64100717
## [19] 7.72669040 7.80652691 7.88052400 7.94868893 8.01102896 8.06755137
## [25] 8.11826344 8.16317242 8.20228560 8.23561024 8.26315362 8.28492300
## [31] 8.30092566 8.31116887 8.31565989 8.31440601 8.30741448 8.29469258
## [37] 8.27624759 8.25208677 8.22221739 8.18664672 8.14538204 8.09843062
## [43] 8.04579973 7.98749663 7.92352860 7.85390291 7.77862683 7.69770763
## [49] 7.61115259 7.51896896 7.42116404 7.31774507 7.20871935 7.09409413
## [55] 6.97387668 6.84807429 6.71669421 6.57974373 6.43723011 6.28916062
## [61] 6.13554253 5.97638312 5.81168965 5.64146939 5.46572963 5.28447762
## [67] 5.09772064 4.90546596 4.70772084 4.50449257 4.29578842 4.08161564
## [73] 3.86198152 3.63689332 3.40635832 3.17038378 2.92897699 2.68214520
## [79] 2.42989569 2.17223573 1.90917259 1.64071354 1.36686585 1.08763680
## [85] 0.80303366 0.51306368 0.21773416 -0.08294765 -0.38897447 -0.70033903
## [91] -1.01703406 -1.33905229 -1.66638645 -1.99902927 -2.33697348 -2.68021179
## [97] -3.02873696 -3.38254170 -3.74161874 -4.10596081summary(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- Perbandingan Model
3.1 Uji ANOVA Bertingkat
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 ' ' 13.2 Perbandingan AIC
AIC(model_linier, model_poli2, model_poli3)## df AIC
## model_linier 3 492.9204
## model_poli2 4 409.4469
## model_poli3 5 411.4307- Visualisasi Perbandingan Model
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()
- Evaluasi dengan RMSE
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- Pemilihan Derajat Optimal via Cross-Validation
Skema 10-fold cross-validation diulang 5 kali untuk menguji derajat polinomial 1 sampai 6.
names(data)## [1] "x" "y" "y_hat_linier" "y_hat_poli2" "y_hat_poli3"
## [6] "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.880719derajat_optimal <- hasil_cv$derajat[which.min(hasil_cv$RMSE)]
cat("Derajat optimal berdasarkan CV:", derajat_optimal, "\n")## Derajat optimal berdasarkan CV: 2##Derajat opimal berdasarkan CV: 2
ggplot(hasil_cv, aes(x = derajat, y = RMSE)) +
geom_line(
color = "yellow",
linewidth = 1
) +
geom_point(
size = 3,
color = "orange"
) +
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()
- Model Final
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-16Catatan
poly(x, derajat, raw = TRUE)digunakan agar koefisien dapat diinterpretasi langsung sebagai \(\beta_1 x + \beta_2 x^2 + \dots\)- Untuk derajat tinggi (\(\geq 4\)),
sebaiknya
xdistandardisasi terlebih dahulu (scale(x)) agar model lebih stabil secara numerik. - Pertimbangkan aturan one-standard-error rule: pilih derajat terkecil dengan RMSE \(\leq\) (RMSE minimum + 1 SD) untuk model yang lebih sederhana namun performanya setara.