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
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)
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
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
AIC(model_linier, model_poli2, model_poli3)
## df AIC
## model_linier 3 492.9204
## model_poli2 4 409.4469
## model_poli3 5 411.4307
library(ggplot2)
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()
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
names(data)
## [1] "x" "y" "pred_linier" "pred_poli2" "pred_poli3"
library(caret)
## Loading required package: lattice
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
## 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()
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