Data simulasi

Data terdiri dari 120 pengamatan. Nilai y dibentuk dari fungsi bergelombang dengan galat bermean 0 dan simpangan baku 0,18.

Ganti nim_seed dengan enam digit angka terakhir NIM sesuai arahan di kelas. Pengaturan seed cukup diubah pada satu baris berikut.

nim_seed <- 123
set.seed(nim_seed)

n <- 120
x <- sort(runif(n, min = 0, max = 1))

f_asli <- function(x) {
  sin(2*pi*x) + 0.35*cos(6*pi*x) + 0.5*x
}

y <- f_asli(x) + rnorm(n, mean = 0, sd = 0.18)
data_prak <- data.frame(x = x, y = y)
grid_pred <- data.frame(x = seq(0, 1, length.out = 500))

Amati pola data

Titik menunjukkan pengamatan. Garis putus-putus menunjukkan fungsi pembangkit yang diketahui karena datanya berupa simulasi.

Diskusi: pada bagian mana pola berubah arah? Apakah kenaikan x selalu diikuti kenaikan y?

Pada data nyata, fungsi pembangkit umumnya tidak diketahui. Pola asosiasi juga belum membuktikan hubungan sebab akibat.

Membentuk kurva

Polinomial menambahkan pangkat x agar model dapat membentuk lengkungan.

Contoh kuadratik ditulis dengan I(x^2). Untuk mencoba derajat lain, gunakan poly().

model_kuadrat <- lm(y ~ x + I(x^2), data = data_prak)

derajat <- 5
model_poly <- lm(y ~ poly(x, degree = derajat), data = data_prak)
pred_poly <- predict(model_poly, newdata = grid_pred)

Model tetap memakai lm() karena linear terhadap koefisien β, meskipun kurvanya melengkung terhadap x.

Secara bawaan, poly() menggunakan basis ortogonal. Koefisien basis ini tidak dibaca langsung sebagai koefisien x, x², dan seterusnya.

Pengaruh derajat

Bandingkan model derajat 5 di atas dengan derajat 1 dan 15 berikut.

derajat_coba <- c(1, 15)
op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), cex = 1)

for (d in derajat_coba) {
  mod <- lm(y ~ poly(x, degree = d), data = data_prak)
  plot_dasar(paste("Derajat", d))
  lines(grid_pred$x, predict(mod, newdata = grid_pred),
        col = "#286CAA", lwd = 2)
}

par(op)

Diskusi: apa yang berubah saat derajat dinaikkan? Perhatikan kelenturan kurva dan perilaku pada bagian tepi.

Kurva yang lebih rumit perlu dinilai dengan validasi sebelum disimpulkan lebih baik.

Prediksi tetap di setiap interval

Rentang x dibagi menjadi beberapa interval. Tanpa peubah lain, prediksi pada setiap interval adalah rerata y dalam interval tersebut.

n_bin <- 8
batas_bin <- seq(0, 1, length.out = n_bin + 1)

data_prak$kelompok <- cut(data_prak$x, breaks = batas_bin,
                          include.lowest = TRUE)
grid_pred$kelompok <- cut(grid_pred$x, breaks = batas_bin,
                          include.lowest = TRUE)

model_step <- lm(y ~ kelompok, data = data_prak)

Dengan K interval, gunakan intersep dan K − 1 indikator. Memasukkan seluruh indikator bersama intersep dapat menghasilkan koefisien NA karena informasi yang berulang sempurna.

Pengaruh jumlah interval

Bandingkan 8 interval di atas dengan 3 dan 20 interval berikut.

bin_coba <- c(3, 20)
op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), cex = 1)

for (b in bin_coba) {
  batas <- seq(0, 1, length.out = b + 1)
  dat_b <- data_prak
  dat_b$kelompok <- cut(dat_b$x, batas, include.lowest = TRUE)
  grid_b <- data.frame(kelompok = cut(grid_pred$x, batas,
                                      include.lowest = TRUE))
  mod <- lm(y ~ kelompok, data = dat_b)
  plot_dasar(paste(b, "interval"))
  lines(grid_pred$x, pred_tangga(mod, grid_b),
        col = "#28856A", lwd = 2, type = "s")
}

par(op)

Diskusi: apa akibat menambah jumlah interval terhadap bentuk fungsi dan banyaknya pengamatan di setiap interval?

Interval kosong tidak memiliki rerata yang dapat diduga. Pada fungsi bantu, prediksi interval tersebut diberi NA dan tampak sebagai celah pada grafik.

Knot sebagai lokasi sambungan

Spline menyambungkan beberapa potongan polinomial pada lokasi yang disebut knot.

Pada cubic spline, nilai fungsi serta turunan pertama dan kedua dibuat kontinu pada setiap knot.

Untuk satu knot pada 0,5, tambahkan basis berikut:

data_prak$h05 <- pmax(data_prak$x - 0.5, 0)^3
grid_pred$h05 <- pmax(grid_pred$x - 0.5, 0)^3

model_spline1 <- lm(y ~ x + I(x^2) + I(x^3) + h05,
                    data = data_prak)

Membandingkan satu dan dua knot

Kandidat kedua memakai knot pada 1/3 dan 2/3. Basis pada data prediksi harus sama dengan basis pada data latih.

data_prak$h1 <- pmax(data_prak$x - 1/3, 0)^3
data_prak$h2 <- pmax(data_prak$x - 2/3, 0)^3
grid_pred$h1 <- pmax(grid_pred$x - 1/3, 0)^3
grid_pred$h2 <- pmax(grid_pred$x - 2/3, 0)^3

model_spline2 <- lm(y ~ x + I(x^2) + I(x^3) + h1 + h2,
                    data = data_prak)

Diskusi: apakah dua knot langsung menghasilkan prediksi yang lebih baik?

Pada perbandingan ini, jumlah sekaligus lokasi knot berubah. Untuk cubic spline biasa dengan K knot interior, terdapat K + 4 koefisien termasuk intersep jika matriks desain berperingkat penuh.

Kehalusan kurva dan penalti

Smoothing spline menyeimbangkan kedekatan terhadap data dengan penalti kelekukan kurva.

• λ lebih besar: kurva lebih halus dan cenderung mendekati garis lurus.
• λ lebih kecil: kurva lebih fleksibel mengikuti data.

model_ss <- smooth.spline(data_prak$x, data_prak$y,
                          all.knots = TRUE, cv = FALSE)

Jika df dan spar belum ditentukan, cv = FALSE memilih kehalusan melalui GCV. Banyak knot tidak membuat kurva wajib melewati seluruh titik.

Mencoba derajat bebas efektif

Derajat bebas efektif (df) yang lebih besar memberi kurva yang lebih fleksibel.

df_coba <- c(4, 30)
op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), cex = 1)

for (d in df_coba) {
  mod <- smooth.spline(data_prak$x, data_prak$y,
                       df = d, all.knots = TRUE)
  plot_dasar(paste("df efektif =", d))
  lines(predict(mod, x = grid_pred$x), col = "#BC4779", lwd = 2)
}

par(op)

Diskusi: mengapa kurva yang sangat fleksibel belum tentu paling baik untuk data baru?

Prinsip 5-fold cross-validation

1. Bagi data menjadi lima bagian.
2. Latih model pada empat bagian dan prediksi satu bagian yang disisihkan.
3. Ulangi sampai setiap bagian pernah menjadi validasi.
4. Bandingkan rata-rata MSE validasi.

Semua kandidat memakai data dan pembagian fold yang sama. Setiap model diestimasi ulang pada data latih di masing-masing fold.

set.seed(nim_seed + 1)
fold_id <- sample(rep(1:5, length.out = nrow(data_prak)))

nama_model <- c("linear", "polinomial", "tangga",
                "spline1", "spline2", "smoothing")
mse_fold <- matrix(NA_real_, nrow = 5, ncol = length(nama_model),
                   dimnames = list(paste0("Fold ", 1:5), nama_model))

for (i in 1:5) {
  latih <- data_prak[fold_id != i, ]
  validasi <- data_prak[fold_id == i, ]

  for (j in seq_along(nama_model)) {
    mod <- fit_model(nama_model[j], latih)
    pred <- pred_model(mod, validasi)
    if (any(!is.finite(pred))) {
      stop("Prediksi validasi belum lengkap. Kurangi jumlah interval, lalu ulangi.")
    }
    mse_fold[i, j] <- mean((validasi$y - pred)^2)
  }
}

round(mse_fold, 4)

model_penuh <- lapply(nama_model, fit_model, dat = data_prak)
names(model_penuh) <- nama_model

mse_latih <- sapply(model_penuh, function(mod) {
  mean((data_prak$y - pred_model(mod, data_prak))^2)
})

hasil_cv <- data.frame(
  Model = nama_model,
  MSE_Latih = as.numeric(mse_latih),
  MSE_CV = colMeans(mse_fold)
)
hasil_cv <- hasil_cv[order(hasil_cv$MSE_CV), ]
rownames(hasil_cv) <- NULL

tabel_tampil <- hasil_cv
tabel_tampil[c("MSE_Latih", "MSE_CV")] <-
  round(tabel_tampil[c("MSE_Latih", "MSE_CV")], 4)
print(tabel_tampil, row.names = FALSE)

Hasil perbandingan

MSE latih memakai pengamatan yang ikut melatih model. MSE CV memakai prediksi pada pengamatan yang disisihkan.

Pada contoh seed 123, MSE CV terkecil diperoleh Cubic spline: 2 knot, sekitar 0.0911.

Diskusi: apakah kandidat dengan MSE latih terkecil juga mempunyai MSE CV terkecil?

Hasil ini berlaku untuk kandidat dan pembagian fold yang dicoba. Jika hasil berdekatan, pertimbangkan kesederhanaan dan kestabilannya pada pembagian ulang.

Kehalusan smoothing spline dipilih ulang dengan GCV pada setiap data latih. Penilaian kinerja akhir setelah pemilihan model idealnya memakai data uji terpisah atau evaluasi luar.

Linear dan polinomial

Fungsi tangga dan cubic spline

Regression spline dan smoothing spline

Diskusi: gabungkan informasi bentuk kurva dan MSE CV untuk menyusun alasan pemilihan model.

Berkas analisis

Sertakan kode R yang dapat dijalankan, grafik utama, tabel perbandingan MSE, dan jawaban latihan. Format serta tenggat pengumpulan mengikuti arahan di kelas.

Jika pengaturan kandidat pada bagian 7 diubah, sesuaikan pembentukan fitur dan fit_model() agar tetap konsisten.