Regresi Polinomial dengan R

#Regresi polinomial adalah model regresi dimana variabel responnya memiliki pangkat kuadratik

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) #sequence artinya berurutan
x

##   [1]  0.0000000  0.1010101  0.2020202  0.3030303  0.4040404  0.5050505
##   [7]  0.6060606  0.7070707  0.8080808  0.9090909  1.0101010  1.1111111
##  [13]  1.2121212  1.3131313  1.4141414  1.5151515  1.6161616  1.7171717
##  [19]  1.8181818  1.9191919  2.0202020  2.1212121  2.2222222  2.3232323
##  [25]  2.4242424  2.5252525  2.6262626  2.7272727  2.8282828  2.9292929
##  [31]  3.0303030  3.1313131  3.2323232  3.3333333  3.4343434  3.5353535
##  [37]  3.6363636  3.7373737  3.8383838  3.9393939  4.0404040  4.1414141
##  [43]  4.2424242  4.3434343  4.4444444  4.5454545  4.6464646  4.7474747
##  [49]  4.8484848  4.9494949  5.0505051  5.1515152  5.2525253  5.3535354
##  [55]  5.4545455  5.5555556  5.6565657  5.7575758  5.8585859  5.9595960
##  [61]  6.0606061  6.1616162  6.2626263  6.3636364  6.4646465  6.5656566
##  [67]  6.6666667  6.7676768  6.8686869  6.9696970  7.0707071  7.1717172
##  [73]  7.2727273  7.3737374  7.4747475  7.5757576  7.6767677  7.7777778
##  [79]  7.8787879  7.9797980  8.0808081  8.1818182  8.2828283  8.3838384
##  [85]  8.4848485  8.5858586  8.6868687  8.7878788  8.8888889  8.9898990
##  [91]  9.0909091  9.1919192  9.2929293  9.3939394  9.4949495  9.5959596
##  [97]  9.6969697  9.7979798  9.8989899 10.0000000

y <- 5 + 2*x - 0.3*x^2 + rnorm(n, mean = 0, sd = 2)
y

##   [1]  3.87904871  4.73860431  8.50921338  5.71952918  6.01768168  9.36370818
##   [7]  7.02376079  3.73403425  5.04655753  5.67892399  9.16227440  7.57147951
##  [13]  7.78501398  7.33033390  6.11666178 10.91542407  8.44443068  3.61650551
##  [19]  9.04733990  6.78781174  5.68039177  7.45661215  5.91095407  6.56945966
##  [25]  6.83532091  4.76404833  9.85892273  8.52988673  5.98053666 10.79198858
##  [31]  9.15871357  7.73094672 10.12052374 10.08960031  9.97343458  9.69837019
##  [37]  9.41362043  8.16053532  7.64488520  7.46219849  6.79393468  7.72160040
##  [43]  5.55460675 12.36515405 10.37888696  5.64634482  7.01026951  6.80008384
##  [49]  9.20455834  7.38300169  7.95536675  7.28450429  7.14260313  9.84617304
##  [55]  6.53192910  9.88479306  3.61660520  7.73947543  6.66797173  6.69603986
##  [61]  6.86120734  4.92893130  5.09269145  3.54136163  3.24821430  5.80601657
##  [67]  5.89641956  4.90092734  6.42825088  8.46656056  3.16088227 -0.70496157
##  [73]  5.68916300  2.01747237  1.81192277  4.98502704  2.10416075 -0.03402802
##  [79]  2.49759319  1.57866050  1.58330676  2.05155254  0.24276319  1.96980605
##  [85] -0.06907230  0.72019080  1.92890913  0.27807648 -1.57778910  1.03192806
##  [91]  0.37543746 -0.86678123 -1.84423840 -3.94176257 -0.33501589 -4.63333215
##  [97]  0.55923899 -1.13894159 -5.07042123 -7.05284180

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

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

head(data.frame(
  Y_asli = data$y,
  Y_topi = predict(model_final),
  Y_residual = residuals(model_final)
), 10)

##      Y_asli   Y_topi   Y_residual
## 1  3.879049 5.339823 -1.460774162
## 2  4.738604 5.519101 -0.780496476
## 3  8.509213 5.692761  2.816452348
## 4  5.719529 5.860804 -0.141274434
## 5  6.017682 6.023229 -0.005546837
## 6  9.363708 6.180036  3.183672422
## 7  7.023761 6.331225  0.692535462
## 8  3.734034 6.476797 -2.742762972
## 9  5.046558 6.616751 -1.570193912
## 10 5.678924 6.751088 -1.072164008

```

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.