1. 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)

data
##              x           y
## 1    0.0000000  3.87904871
## 2    0.1010101  4.73860431
## 3    0.2020202  8.50921338
## 4    0.3030303  5.71952918
## 5    0.4040404  6.01768168
## 6    0.5050505  9.36370818
## 7    0.6060606  7.02376079
## 8    0.7070707  3.73403425
## 9    0.8080808  5.04655753
## 10   0.9090909  5.67892399
## 11   1.0101010  9.16227440
## 12   1.1111111  7.57147951
## 13   1.2121212  7.78501398
## 14   1.3131313  7.33033390
## 15   1.4141414  6.11666178
## 16   1.5151515 10.91542407
## 17   1.6161616  8.44443068
## 18   1.7171717  3.61650551
## 19   1.8181818  9.04733990
## 20   1.9191919  6.78781174
## 21   2.0202020  5.68039177
## 22   2.1212121  7.45661215
## 23   2.2222222  5.91095407
## 24   2.3232323  6.56945966
## 25   2.4242424  6.83532091
## 26   2.5252525  4.76404833
## 27   2.6262626  9.85892273
## 28   2.7272727  8.52988673
## 29   2.8282828  5.98053666
## 30   2.9292929 10.79198858
## 31   3.0303030  9.15871357
## 32   3.1313131  7.73094672
## 33   3.2323232 10.12052374
## 34   3.3333333 10.08960031
## 35   3.4343434  9.97343458
## 36   3.5353535  9.69837019
## 37   3.6363636  9.41362043
## 38   3.7373737  8.16053532
## 39   3.8383838  7.64488520
## 40   3.9393939  7.46219849
## 41   4.0404040  6.79393468
## 42   4.1414141  7.72160040
## 43   4.2424242  5.55460675
## 44   4.3434343 12.36515405
## 45   4.4444444 10.37888696
## 46   4.5454545  5.64634482
## 47   4.6464646  7.01026951
## 48   4.7474747  6.80008384
## 49   4.8484848  9.20455834
## 50   4.9494949  7.38300169
## 51   5.0505051  7.95536675
## 52   5.1515152  7.28450429
## 53   5.2525253  7.14260313
## 54   5.3535354  9.84617304
## 55   5.4545455  6.53192910
## 56   5.5555556  9.88479306
## 57   5.6565657  3.61660520
## 58   5.7575758  7.73947543
## 59   5.8585859  6.66797173
## 60   5.9595960  6.69603986
## 61   6.0606061  6.86120734
## 62   6.1616162  4.92893130
## 63   6.2626263  5.09269145
## 64   6.3636364  3.54136163
## 65   6.4646465  3.24821430
## 66   6.5656566  5.80601657
## 67   6.6666667  5.89641956
## 68   6.7676768  4.90092734
## 69   6.8686869  6.42825088
## 70   6.9696970  8.46656056
## 71   7.0707071  3.16088227
## 72   7.1717172 -0.70496157
## 73   7.2727273  5.68916300
## 74   7.3737374  2.01747237
## 75   7.4747475  1.81192277
## 76   7.5757576  4.98502704
## 77   7.6767677  2.10416075
## 78   7.7777778 -0.03402802
## 79   7.8787879  2.49759319
## 80   7.9797980  1.57866050
## 81   8.0808081  1.58330676
## 82   8.1818182  2.05155254
## 83   8.2828283  0.24276319
## 84   8.3838384  1.96980605
## 85   8.4848485 -0.06907230
## 86   8.5858586  0.72019080
## 87   8.6868687  1.92890913
## 88   8.7878788  0.27807648
## 89   8.8888889 -1.57778910
## 90   8.9898990  1.03192806
## 91   9.0909091  0.37543746
## 92   9.1919192 -0.86678123
## 93   9.2929293 -1.84423840
## 94   9.3939394 -3.94176257
## 95   9.4949495 -0.33501589
## 96   9.5959596 -4.63333215
## 97   9.6969697  0.55923899
## 98   9.7979798 -1.13894159
## 99   9.8989899 -5.07042123
## 100 10.0000000 -7.05284180

2. 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)

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

3. 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 ' ' 1

3.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

4. 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()

5. 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

6. 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"           "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()

7. 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-16

Catatan

  • 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 x distandardisasi 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.