Modul Praktikum

Regresi Nonlinear

1 Pengantar

Pada praktikum sebelumnya, kita menggunakan garis lurus untuk menjelaskan hubungan antara peubah penjelas dan peubah respons. Namun, pola hubungan pada data tidak selalu lurus. Ada hubungan yang melengkung, ada pula yang berubah lebih cepat pada rentang tertentu.

Kali ini kita akan mencoba empat cara untuk menggambarkan pola tersebut: regresi polinomial, fungsi tangga, regression spline, dan smoothing spline. Kita akan menggunakan data simulasi agar bentuk hubungan yang sebenarnya diketahui. Dengan begitu, kita bisa melihat bagian mana yang berhasil ditangkap oleh setiap model.

Setelah mengikuti praktikum ini, mahasiswa diharapkan dapat:

1. mengenali pola hubungan yang kurang tepat dijelaskan dengan garis lurus;
2. membentuk keempat model di R;
3. membaca hasil pendugaan dan menghubungkannya dengan grafik;
4. membandingkan kemampuan prediksi model menggunakan cross-validation.

install.packages(c("rmarkdown", "knitr"))

2 Mengenali hubungan yang melengkung

Bayangkan kita ingin menjelaskan hubungan antara biaya iklan dan penjualan. Pada awalnya, tambahan biaya iklan mungkin diikuti kenaikan penjualan yang cukup besar. Setelah biaya iklan semakin tinggi, tambahan penjualannya bisa semakin kecil. Hubungan seperti ini sulit digambarkan dengan satu garis lurus karena laju perubahannya tidak tetap.

Pada regresi linear sederhana, model yang digunakan adalah

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i. \]

Koefisien \(\beta_1\) menyatakan perubahan rata-rata respons untuk setiap kenaikan satu satuan \(x\). Besarnya sama di sepanjang rentang \(x\). Jika hubungan pada data melengkung, kita perlu memberi ruang agar laju perubahan tersebut dapat berbeda.

Bentuk umumnya dapat ditulis sebagai

\[ y_i = f(x_i) + \varepsilon_i, \]

dengan \(f(x_i)\) menyatakan pola hubungan yang ingin kita duga, sedangkan \(\varepsilon_i\) menyatakan galat.

Metode	Cara membentuk hubungan
Regresi polinomial	Menambahkan pangkat \(x\), misalnya \(x^2\) dan \(x^3\).
Fungsi tangga	Membagi \(x\) menjadi interval, lalu menduga satu nilai rata-rata pada setiap interval.
Regression spline	Menggunakan polinomial per bagian yang disambungkan secara mulus pada titik tertentu.
Smoothing spline	Membentuk kurva dengan menyeimbangkan kedekatan terhadap data dan kehalusan kurva.

Istilah nonlinear pada modul ini merujuk pada bentuk hubungan antara \(x\) dan rata-rata \(y\). Regresi polinomial dan regression spline yang kita gunakan tetap linear terhadap koefisiennya, sehingga dapat diduga dengan fungsi lm().

Model regresi menunjukkan pola hubungan dalam data. Hubungan tersebut belum cukup untuk menyimpulkan bahwa perubahan \(x\) menyebabkan perubahan \(y\).

3 Menyiapkan data dan grafik

3.1 Membuat data simulasi

Kita mulai dengan hubungan kuadratik:

\[ y_i = 5 + 2x_i + 3x_i^2 + \varepsilon_i, \qquad \varepsilon_i \sim N(0,1). \]

Nilai 5, 2, dan 3 merupakan koefisien yang kita tentukan sendiri. Galat ditambahkan agar titik-titik data tidak seluruhnya berada tepat pada kurva.

set.seed(123)

data.x <- rnorm(1000, mean = 1, sd = 1)
err <- rnorm(1000, mean = 0, sd = 1)
y <- 5 + 2 * data.x + 3 * data.x^2 + err

dat_awal <- data.frame(data.x, y)
summary(dat_awal)

##      data.x             y        
##  Min.   :-1.810   Min.   : 2.73  
##  1st Qu.: 0.372   1st Qu.: 6.48  
##  Median : 1.009   Median :10.19  
##  Mean   : 1.016   Mean   :13.12  
##  3rd Qu.: 1.665   3rd Qu.:16.86  
##  Max.   : 4.241   Max.   :67.33

set.seed(123) membuat hasil simulasi dapat diulang. Sementara itu, rnorm() digunakan untuk membangkitkan bilangan acak dari sebaran normal. Pada kode di atas, kita membangkitkan 1.000 nilai \(x\) serta 1.000 galat.

3.2 Menyiapkan fungsi untuk menggambar

Kita akan menggambar beberapa grafik dengan tampilan yang sama. Agar tidak menulis pengaturan grafik berulang kali, kita simpan pengaturan tersebut dalam fungsi berikut.

pal <- c(
  linear = "#D64550", polinomial = "#2F6690",
  tangga = "#3A7D44", spline = "#D48600",
  smoothing = "#7B4B94", asli = "#202020"
)

plot_sebaran <- function(x, y, judul = "") {
  plot(
    x, y, type = "n",
    xlab = "x", ylab = "y", main = judul
  )
  grid(col = "grey90", lty = 1)
  points(x, y, pch = 16,
         col = adjustcolor("grey30", alpha.f = 0.3))
}

plot_resid <- function(model, judul, warna) {
  plot(
    fitted(model), resid(model),
    pch = 16, col = adjustcolor("grey30", alpha.f = 0.3),
    xlab = "Nilai prediksi", ylab = "Residual", main = judul
  )
  abline(h = 0, lty = 2)
  lines(lowess(fitted(model), resid(model)),
        col = warna, lwd = 2)
}

Fungsi pertama membuat diagram pencar. Fungsi kedua menggambarkan residual terhadap nilai prediksi. fitted() mengambil nilai dugaan pada data yang digunakan untuk menduga model, sedangkan resid() mengambil residualnya.

plot_sebaran(dat_awal$data.x, dat_awal$y, "Pola hubungan pada data simulasi")

Perhatikan pola titik-titiknya. Pada bagian kanan, nilai \(y\) semakin cepat meningkat ketika \(x\) bertambah. Satu garis lurus kemungkinan belum cukup untuk mengikuti pola ini.

Untuk menggambar garis model, kita juga menyiapkan urutan nilai \(x\) dari kecil ke besar.

ix_awal <- order(dat_awal$data.x)

Pengurutan diperlukan karena lines() menyambungkan titik sesuai urutan data yang diberikan. Tanpa pengurutan, garis dapat bergerak bolak-balik mengikuti urutan data acak.

4 Regresi linear sebagai pembanding

4.1 Pengertian regresi linear

Regresi linear sederhana adalah metode untuk menjelaskan hubungan antara satu peubah penjelas dan rata-rata peubah respons menggunakan garis lurus. Dalam praktikum ini, peubah penjelas diberi simbol \(x\), sedangkan respons diberi simbol \(y\).

Ciri utamanya adalah kemiringan garis yang tetap. Menurut model, kenaikan satu satuan \(x\) diikuti perubahan rata-rata \(y\) dengan besar yang sama pada setiap posisi \(x\). Jika kemiringannya positif, garis naik; jika negatif, garis turun.

Sebagai contoh, model dapat menggambarkan hubungan jam belajar dan nilai ujian dengan garis lurus ketika tambahan rata-rata nilai untuk setiap jam belajar dianggap tetap pada rentang yang diamati. Kesesuaian anggapan ini tetap perlu diperiksa melalui data.

Pada praktikum ini, regresi linear digunakan sebagai model awal. Kita akan melihat apakah garis lurus sudah cukup mengikuti pola data sebelum mencoba bentuk yang lebih lentur.

4.2 Menduga model

Kita pasangkan regresi linear terlebih dahulu. Model ini menjadi pembanding untuk melihat apakah penambahan kelengkungan memang membantu menjelaskan data.

lin.mod <- lm(y ~ data.x, data = dat_awal)
summary(lin.mod)

## 
## Call:
## lm(formula = y ~ data.x, data = dat_awal)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -5.69  -2.57  -1.43   1.20  27.18 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.606      0.190    24.2   <2e-16 ***
## data.x         8.379      0.134    62.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.2 on 998 degrees of freedom
## Multiple R-squared:  0.797,  Adjusted R-squared:  0.796 
## F-statistic: 3.91e+03 on 1 and 998 DF,  p-value: <2e-16

Rumus y ~ data.x berarti kita menjelaskan \(y\) menggunakan \(x\). R secara otomatis memasukkan intersep.

Pada tabel koefisien, Estimate berisi nilai dugaan koefisien, Std. Error menunjukkan galat bakunya, dan Pr(>|t|) berisi p-value uji koefisien. Bagian bawah keluaran memuat simpangan baku residual serta \(R^2\).

4.3 Melihat garis hasil pendugaan

plot_sebaran(dat_awal$data.x, dat_awal$y, "Hasil regresi linear")
lines(
  dat_awal$data.x[ix_awal], fitted(lin.mod)[ix_awal],
  col = pal["linear"], lwd = 3
)

Persamaan garis yang diperoleh adalah

\[ \widehat{y} = \text{4.61} + \text{8.38}\,x. \]

Menurut model ini, setiap kenaikan satu satuan \(x\) diikuti kenaikan rata-rata \(y\) sebesar sekitar 8.38 satuan. Ini merupakan ringkasan berdasarkan garis lurus yang dipasangkan; kita masih perlu memeriksa apakah bentuk tersebut sesuai dengan data.

Nilai \(R^2\) model adalah 0.797. Meskipun angkanya cukup tinggi, garis belum mengikuti lengkungan titik-titik dengan baik. Pada beberapa bagian garis berada di atas pola data, sedangkan pada bagian lainnya berada di bawahnya.

Nilai \(R^2\) yang tinggi belum menjamin bentuk hubungan sudah tepat. Selalu hubungkan keluaran model dengan grafiknya.

5 Regresi polinomial

5.1 Pengertian regresi polinomial

Regresi polinomial adalah metode regresi yang menggunakan pangkat-pangkat peubah penjelas untuk menggambarkan hubungan yang melengkung antara \(x\) dan rata-rata \(y\). Selain suku \(x\), model dapat memuat \(x^2\), \(x^3\), dan seterusnya.

Pangkat tertinggi yang digunakan disebut derajat polinomial. Model derajat dua menggunakan suku sampai \(x^2\) dan dapat membentuk kurva seperti huruf U atau U terbalik. Model derajat tiga menggunakan suku sampai \(x^3\), sehingga dapat mengikuti bentuk lengkungan yang lebih beragam.

Sebagai contoh, hubungan dosis pupuk dan rata-rata hasil panen mungkin meningkat pada dosis awal, lalu kenaikannya melambat atau hasilnya menurun pada dosis yang berlebihan. Pada rentang tertentu, pola tersebut dapat didekati dengan kurva polinomial.

Satu persamaan polinomial digunakan untuk seluruh rentang \(x\). Derajat yang lebih tinggi memberi lebih banyak kelenturan, tetapi belum tentu membuat prediksi lebih baik. Di sini kita mulai dari derajat dua karena data simulasi memiliki pola kuadratik.

5.2 Menambahkan suku kuadrat

Regresi polinomial menambahkan pangkat dari \(x\). Untuk polinomial derajat dua, modelnya adalah

\[ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \varepsilon_i. \]

Tambahan \(x^2\) memungkinkan kurva melengkung. Model ini masih linear terhadap \(\beta_0,\beta_1,\beta_2\), sehingga kita tetap menggunakan lm().

pol.mod <- lm(y ~ data.x + I(data.x^2), data = dat_awal)
summary(pol.mod)

## 
## Call:
## lm(formula = y ~ data.x + I(data.x^2), data = dat_awal)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.032 -0.694  0.005  0.712  3.286 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.9519     0.0457     108   <2e-16 ***
## data.x        2.1073     0.0586      36   <2e-16 ***
## I(data.x^2)   2.9908     0.0234     128   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.01 on 997 degrees of freedom
## Multiple R-squared:  0.988,  Adjusted R-squared:  0.988 
## F-statistic: 4.22e+04 on 2 and 997 DF,  p-value: <2e-16

Penulisan I(data.x^2) memberi tahu R agar menghitung kuadrat \(x\) sebagai operasi aritmetika biasa. Ini diperlukan karena simbol di dalam rumus model memiliki aturan khusus.

5.3 Membandingkan kurva dan residual

plot_sebaran(dat_awal$data.x, dat_awal$y,
             "Regresi linear dan polinomial derajat dua")
lines(dat_awal$data.x[ix_awal], fitted(lin.mod)[ix_awal],
      col = pal["linear"], lwd = 3)
lines(dat_awal$data.x[ix_awal], fitted(pol.mod)[ix_awal],
      col = pal["polinomial"], lwd = 3)
legend(
  "topleft", c("Linear", "Polinomial derajat dua"),
  col = pal[c("linear", "polinomial")], lwd = 3, bty = "n"
)

Kurva polinomial lebih mampu mengikuti pola titik-titik. Untuk memeriksanya lebih lanjut, kita lihat residual, yaitu selisih antara respons yang diamati dan nilai dugaannya:

\[ e_i = y_i - \widehat{y}_i. \]

par(mfrow = c(1, 2))
plot_resid(lin.mod, "Residual model linear", pal["linear"])
plot_resid(pol.mod, "Residual model polinomial", pal["polinomial"])

par(mfrow = c(1, 1))

Garis berwarna pada grafik merupakan garis penghalus untuk membantu melihat pola residual. Pada model linear, residual masih memperlihatkan pola sistematis. Pada model polinomial, polanya berkurang dan residual lebih menyebar di sekitar nol.

Grafik ini membantu memeriksa kecocokan bentuk rata-rata dan perubahan sebaran residual. Pemeriksaan asumsi lain, seperti independensi dan normalitas galat untuk keperluan inferensi, memerlukan pemeriksaan tambahan.

5.4 Apakah penambahan suku kuadrat membantu?

Kedua model bersarang: model linear merupakan model polinomial ketika \(\beta_2=0\). Karena menggunakan respons dan pengamatan yang sama, kita dapat membandingkannya dengan uji F berikut.

anova(lin.mod, pol.mod)

## Analysis of Variance Table
## 
## Model 1: y ~ data.x
## Model 2: y ~ data.x + I(data.x^2)
##   Res.Df   RSS Df Sum of Sq     F Pr(>F)    
## 1    998 17601                              
## 2    997  1011  1     16590 16366 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hipotesis nolnya adalah \(H_0:\beta_2=0\). P-value yang kecil memberikan bukti bahwa penambahan suku kuadrat memperbaiki kecocokan terhadap data, dengan asumsi uji F terpenuhi. Hasil ini belum menggantikan evaluasi prediksi pada data yang tidak digunakan untuk menduga model.

Bandingkan juga koefisien hasil pendugaan dengan nilai yang digunakan saat simulasi.

knitr::kable(
  data.frame(
    Koefisien = c("Intersep", "x", "x kuadrat"),
    Nilai_simulasi = c(5, 2, 3),
    Hasil_pendugaan = unname(coef(pol.mod))
  ),
  digits = 3, row.names = FALSE
)

Koefisien	Nilai_simulasi	Hasil_pendugaan
Intersep	5	4.952
x	2	2.107
x kuadrat	3	2.991

Nilai dugaannya tidak harus sama persis dengan koefisien simulasi karena data mengandung galat acak.

5.5 Membaca koefisien polinomial

Pada model linear sederhana, satu koefisien menyatakan kemiringan garis di seluruh rentang \(x\). Pada model kuadratik, kemiringannya bergantung pada posisi \(x\):

\[ \frac{d\,\widehat{y}}{dx} = \widehat{\beta}_1 + 2\widehat{\beta}_2 x. \]

Artinya, koefisien \(x\) tidak dibaca sendirian sebagai perubahan yang tetap untuk setiap kenaikan satu satuan \(x\). Ketika \(x\) berubah, \(x^2\) juga berubah.

b <- coef(pol.mod)
x0 <- c(0, 1, 2)

knitr::kable(
  data.frame(
    x = x0,
    Kemiringan_kurva =
      b[["data.x"]] + 2 * b[["I(data.x^2)"]] * x0
  ),
  digits = 2, row.names = FALSE
)

x	Kemiringan_kurva
0	2.11
1	8.09
2	14.07

Tabel tersebut menunjukkan laju perubahan sesaat di setiap posisi \(x\). Untuk perubahan tepat dari \(x\) ke \(x+1\), selisih nilai prediksinya adalah \(\widehat{\beta}_1+\widehat{\beta}_2(2x+1)\), sehingga berbeda dari turunan pada satu titik.

6 Regresi fungsi tangga

6.1 Pengertian regresi fungsi tangga

Regresi fungsi tangga adalah metode regresi yang membagi rentang peubah penjelas menjadi beberapa interval, lalu menggunakan satu nilai rata-rata respons sebagai prediksi untuk setiap interval. Karena nilai prediksinya tetap di dalam interval dan dapat berubah ketika melewati batas interval, grafiknya menyerupai anak tangga.

Cara kerjanya adalah mengelompokkan nilai \(x\) terlebih dahulu. Semua pengamatan dalam kelompok yang sama memperoleh prediksi yang sama. Jadi, model menangkap perbedaan rata-rata antarkelompok, tetapi tidak menggambarkan perubahan respons secara bertahap di dalam kelompok.

Sebagai contoh, usia dapat dibagi menjadi kelompok 18–25 tahun, 26–35 tahun, dan 36–45 tahun. Jika usia menjadi satu-satunya penjelas, model memberikan satu prediksi rata-rata pengeluaran untuk setiap kelompok usia tersebut. Dua orang yang berbeda usia tetapi masih berada dalam kelompok yang sama mendapatkan prediksi yang sama.

Titik yang memisahkan interval disebut batas interval atau titik potong. Jumlah dan letak batas ini memengaruhi hasil model. Fungsi tangga berguna ketika perbedaan antarkelompok memang menjadi perhatian, tetapi perubahan mendadak pada batas kelompok belum tentu mencerminkan hubungan yang sebenarnya.

6.2 Membagi nilai x menjadi interval

Sekarang kita mencoba pendekatan lain. Rentang \(x\) dibagi menjadi tiga interval:

1. \(x \leq 0\);
2. \(0 < x \leq 2\);
3. \(x > 2\).

Setiap interval memiliki satu nilai rata-rata dugaan. Dengan demikian, dua pengamatan yang berada dalam interval sama akan mendapatkan prediksi yang sama, meskipun nilai \(x\)-nya berbeda.

dat_awal$kelas <- cut(
  dat_awal$data.x,
  breaks = c(-Inf, 0, 2, Inf),
  labels = c("x <= 0", "0 < x <= 2", "x > 2")
)

table(dat_awal$kelas)

## 
##     x <= 0 0 < x <= 2      x > 2 
##        165        678        157

Fungsi cut() mengubah nilai numerik menjadi kategori berdasarkan batas interval. Dengan pengaturan bawaan, batas kanan masuk ke dalam interval: nilai 0 masuk ke kelompok pertama dan nilai 2 masuk ke kelompok kedua.

6.3 Menduga dan membaca model

step.mod <- lm(y ~ kelas, data = dat_awal)
summary(step.mod)

## 
## Call:
## lm(formula = y ~ kelas, data = dat_awal)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -10.39  -3.53  -0.53   2.53  36.88 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        5.210      0.399    13.1   <2e-16 ***
## kelas0 < x <= 2    5.821      0.445    13.1   <2e-16 ***
## kelasx > 2        25.239      0.571    44.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.12 on 997 degrees of freedom
## Multiple R-squared:  0.698,  Adjusted R-squared:  0.697 
## F-statistic: 1.15e+03 on 2 and 997 DF,  p-value: <2e-16

Kelompok pertama, \(x\leq0\), digunakan sebagai kategori pembanding. Modelnya dapat ditulis sebagai

\[ y_i = \beta_0 + \beta_1 I(0<x_i\leq2) + \beta_2 I(x_i>2) + \varepsilon_i. \]

Fungsi indikator \(I(\cdot)\) bernilai 1 jika syarat di dalamnya terpenuhi dan 0 jika tidak.

Cara membaca koefisiennya adalah sebagai berikut:

• intersep menyatakan rata-rata \(y\) pada kelompok pertama;
• koefisien kelompok kedua menyatakan selisih rata-ratanya terhadap kelompok pertama;
• koefisien kelompok ketiga menyatakan selisih rata-ratanya terhadap kelompok pertama.

Nilai rata-rata dugaan untuk kelompok kedua dan ketiga diperoleh dengan menjumlahkan intersep dan koefisien kelompok tersebut.

rata_kelompok <- aggregate(y ~ kelas, data = dat_awal, FUN = mean)
names(rata_kelompok)[2] <- "Rata_rata_y"
knitr::kable(rata_kelompok, digits = 3, row.names = FALSE)

kelas	Rata_rata_y
x <= 0	5.21
0 < x <= 2	11.03
x > 2	30.45

Pada regresi yang hanya menggunakan kategori ini, nilai prediksi untuk setiap kelompok sama dengan rata-rata respons dalam kelompoknya.

6.4 Menggambar fungsi tangga

Untuk menunjukkan bahwa prediksi konstan di setiap interval, kita menggambar tiga garis mendatar secara terpisah.

plot_sebaran(dat_awal$data.x, dat_awal$y, "Fungsi tangga dengan tiga interval")
abline(v = c(0, 2), col = "grey55", lty = 2)

batas <- c(min(dat_awal$data.x), 0, 2, max(dat_awal$data.x))
rerata <- tapply(dat_awal$y, dat_awal$kelas, mean)

segments(
  x0 = batas[1:3], y0 = as.numeric(rerata),
  x1 = batas[2:4], y1 = as.numeric(rerata),
  col = pal["tangga"], lwd = 3
)

Fungsi tangga mudah dijelaskan, tetapi bentuknya berubah secara mendadak pada batas interval. Di dalam interval, perubahan \(x\) tidak diikuti perubahan prediksi. Pada data yang hubungannya terus melengkung, pembagian menjadi tiga interval dapat terlalu kasar.

Hasil model juga bergantung pada batas yang dipilih. Jika tujuannya memang membandingkan kelompok, pendekatan ini dapat membantu. Jika tujuannya menangkap perubahan yang mulus, kita bisa mempertimbangkan spline.

6.5 Jika variabel dummy dibuat sendiri

Satu kategori dengan tiga kelompok cukup diwakili dua dummy jika model memiliki intersep. Contohnya:

dat_awal$c2 <- as.integer(dat_awal$data.x > 0 & dat_awal$data.x <= 2)
dat_awal$c3 <- as.integer(dat_awal$data.x > 2)

step.manual <- lm(y ~ c2 + c3, data = dat_awal)

all.equal(
  unname(fitted(step.mod)),
  unname(fitted(step.manual))
)

## [1] TRUE

Kelompok pertama diwakili kondisi ketika kedua dummy bernilai nol. Jika tiga dummy dimasukkan sekaligus bersama intersep, kolom-kolom tersebut saling bergantung dan model tidak dapat menduga seluruh koefisien secara unik. R biasanya menandai salah satu koefisien sebagai NA. Penggunaan cut() sebagai faktor memudahkan kita menghindari masalah ini.

7 Regression spline

7.1 Pengertian regression spline

Regression spline adalah metode regresi yang membagi rentang peubah penjelas menjadi beberapa bagian, memasangkan fungsi polinomial pada bagian-bagian tersebut, lalu menyambungkannya dengan syarat tertentu. Titik tempat bagian-bagian kurva bertemu disebut knot atau titik simpul.

Berbeda dari satu polinomial yang berlaku untuk seluruh rentang data, spline memberi ruang agar bentuk hubungan dapat berbeda pada bagian yang berbeda. Berbeda pula dari fungsi tangga, prediksi spline dapat berubah di dalam setiap bagian. Untuk spline kubik yang digunakan di sini, sambungannya dibuat mulus: nilai fungsi serta turunan pertama dan keduanya kontinu di knot.

Sebagai contoh, rata-rata tinggi tanaman dapat bertambah lambat pada awal pengamatan, lebih cepat pada masa pertumbuhan, lalu kembali melambat. Spline dapat digunakan untuk mendekati pola perubahan tersebut tanpa memaksakan satu bentuk lengkungan sederhana untuk seluruh waktu pengamatan.

Pada regression spline, kita menentukan jumlah dan lokasi knot. Pilihan ini mengatur kelenturan model. Dalam praktikum ini kita menggunakan polinomial derajat tiga, sehingga modelnya disebut spline kubik, lalu membandingkan hasil penggunaan satu dan dua knot.

7.2 Menyiapkan data dengan galat lebih besar

Pada bagian ini, simpangan baku galat kita naikkan dari 1 menjadi 5 agar perbedaan kelenturan kurva lebih mudah diamati. Bentuk hubungan dan nilai \(x\) tetap sama.

Data baru disimpan dalam objek dat, sedangkan dat_awal tetap tersedia. Dengan demikian, hasil dari dua simulasi tidak tertukar.

dat <- data.frame(
  data.x = dat_awal$data.x,
  y = 5 + 2 * dat_awal$data.x + 3 * dat_awal$data.x^2 + 5 * err
)
ix <- order(dat$data.x)

plot_sebaran(dat$data.x, dat$y, "Data simulasi dengan simpangan baku galat 5")

Model pembangkit data sekarang adalah

\[ y_i=5+2x_i+3x_i^2+\varepsilon_i, \qquad \varepsilon_i\sim N(0,25). \]

Jangan langsung membandingkan \(R^2\), residual, atau galat model pada data ini dengan model sebelumnya. Responsnya sudah berubah. Pada bagian perbandingan nanti, semua model akan diduga ulang menggunakan dat.

7.3 Mengapa perlu titik sambung?

Regression spline membagi rentang \(x\) menjadi beberapa bagian. Pada setiap bagian, hubungan dijelaskan oleh polinomial. Bagian-bagian tersebut disambungkan pada titik yang disebut knot.

Sebagai contoh, gunakan knot di \(x=1\). Kita coba menduga polinomial kubik di kiri dan kanan knot secara terpisah.

dt1 <- subset(dat, data.x <= 1)
dt2 <- subset(dat, data.x > 1)

cub.mod1 <- lm(y ~ data.x + I(data.x^2) + I(data.x^3), data = dt1)
cub.mod2 <- lm(y ~ data.x + I(data.x^2) + I(data.x^3), data = dt2)

plot_sebaran(dat$data.x, dat$y, "Dua polinomial kubik yang diduga terpisah")
abline(v = 1, col = "grey55", lty = 2)

x_kiri <- seq(min(dat$data.x), 1, length.out = 150)
x_kanan <- seq(1, max(dat$data.x), length.out = 150)

lines(x_kiri, predict(cub.mod1, data.frame(data.x = x_kiri)),
      col = pal["polinomial"], lwd = 3)
lines(x_kanan, predict(cub.mod2, data.frame(data.x = x_kanan)),
      col = pal["spline"], lwd = 3)
legend("topleft", c("Bagian kiri", "Bagian kanan"),
       col = pal[c("polinomial", "spline")], lwd = 3, bty = "n")

Karena kedua model diduga sendiri-sendiri, tidak ada syarat yang memaksa keduanya bertemu di knot. Kita dapat memeriksa selisih prediksinya tepat pada \(x=1\).

kiri <- predict(cub.mod1, newdata = data.frame(data.x = 1))
kanan <- predict(cub.mod2, newdata = data.frame(data.x = 1))

data.frame(
  Prediksi_kiri = unname(kiri),
  Prediksi_kanan = unname(kanan),
  Selisih = unname(kanan - kiri)
)

##   Prediksi_kiri Prediksi_kanan Selisih
## 1         9.634          11.12    1.49

Jika nilai kedua prediksi berbeda, terdapat loncatan pada titik sambung. Untuk hubungan yang diharapkan mulus, kita perlu memberi syarat pada sambungan tersebut.

7.4 Membentuk spline kubik

Pada spline kubik, kurva beserta turunan pertama dan keduanya dibuat kontinu di knot. Salah satu cara membentuknya adalah menggunakan basis pangkat terpotong:

\[ y_i = \beta_0+\beta_1x_i+\beta_2x_i^2+\beta_3x_i^3 +\beta_4(x_i-\xi)_+^3+\varepsilon_i, \]

dengan

\[ (x-\xi)_+^3= \begin{cases} 0, & x\leq\xi,\\ (x-\xi)^3, & x>\xi. \end{cases} \]

Kita tetap memakai \(\xi=1\). Sebelum knot, suku tambahan tersebut bernilai nol. Setelah knot, nilainya mengikuti \((x-1)^3\). Suku ini memungkinkan bentuk kurva berubah sambil mempertahankan sambungan yang mulus.

dat$hx <- ifelse(dat$data.x > 1, (dat$data.x - 1)^3, 0)
head(dat)

##   data.x      y        hx
## 1 0.4395  1.480 0.0000000
## 2 0.7698  3.118 0.0000000
## 3 2.5587 29.668 3.7869935
## 4 1.0705  9.918 0.0003505
## 5 1.1293 -1.662 0.0021611
## 6 2.7151 37.748 5.0447743

cubspline.mod <- lm(
  y ~ data.x + I(data.x^2) + I(data.x^3) + hx,
  data = dat
)
summary(cubspline.mod)

## 
## Call:
## lm(formula = y ~ data.x + I(data.x^2) + I(data.x^3) + hx, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.404  -3.440  -0.017   3.570  16.603 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.168      0.301   17.16  < 2e-16 ***
## data.x         2.170      0.411    5.28  1.6e-07 ***
## I(data.x^2)    2.265      0.355    6.39  2.6e-10 ***
## I(data.x^3)    0.461      0.241    1.92    0.055 .  
## hx            -0.653      0.421   -1.55    0.121    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.03 on 995 degrees of freedom
## Multiple R-squared:  0.784,  Adjusted R-squared:  0.783 
## F-statistic:  902 on 4 and 995 DF,  p-value: <2e-16

plot_sebaran(dat$data.x, dat$y, "Spline kubik dengan knot di x = 1")
abline(v = 1, col = "grey55", lty = 2)
lines(dat$data.x[ix], fitted(cubspline.mod)[ix],
      col = pal["spline"], lwd = 3)

Kurva kini tersambung mulus di \(x=1\). Koefisien hx mengatur tambahan komponen kubik setelah knot. P-value koefisien ini menguji kontribusi suku tambahan tersebut dengan mempertimbangkan suku lain dalam model. Hasil yang tidak signifikan tidak membuktikan bahwa seluruh bentuk hubungan di sekitar knot sama.

Untuk membaca spline, utamakan bentuk kurva dan hasil evaluasi prediksinya. Koefisien basisnya biasanya kurang mudah dijelaskan satu per satu dalam konteks substantif.

7.5 Mencoba dua knot

Kita bandingkan model tadi dengan spline yang memiliki dua knot, yaitu di \(x=0\) dan \(x=2\).

cubspline.mod2 <- lm(
  y ~ data.x + I(data.x^2) + I(data.x^3) +
    I(pmax(data.x, 0)^3) + I(pmax(data.x - 2, 0)^3),
  data = dat
)

par(mfrow = c(1, 2))
plot_sebaran(dat$data.x, dat$y, "Satu knot: x = 1")
abline(v = 1, col = "grey55", lty = 2)
lines(dat$data.x[ix], fitted(cubspline.mod)[ix],
      col = pal["spline"], lwd = 3)

plot_sebaran(dat$data.x, dat$y, "Dua knot: x = 0 dan x = 2")
abline(v = c(0, 2), col = "grey55", lty = 2)
lines(dat$data.x[ix], fitted(cubspline.mod2)[ix],
      col = pal["spline"], lwd = 3)

par(mfrow = c(1, 1))

Penulisan pmax(data.x - 1, 0)^3 menghasilkan basis yang sama dengan ifelse(data.x > 1, (data.x - 1)^3, 0). Jika ditulis langsung di dalam rumus model, basis juga dihitung otomatis ketika kita melakukan prediksi untuk data baru.

Semakin banyak knot, semakin besar ruang bagi kurva untuk berubah. Spline kubik dengan \(K\) knot interior memiliki \(K+4\) koefisien, termasuk intersep, apabila matriks model berpangkat penuh. Knot dapat ditempatkan pada jarak yang sama atau berdasarkan kuantil \(x\). Pemilihan berdasarkan kuantil membantu menyebarkan pengamatan di antara bagian-bagian kurva.

Penambahan knot belum tentu memperbaiki prediksi. Grafik memberi gambaran awal, sedangkan cross-validation akan membantu menilai hasilnya.

8 Smoothing spline

8.1 Pengertian smoothing spline

Smoothing spline adalah metode regresi yang menduga hubungan antara peubah penjelas dan rata-rata respons dengan kurva mulus, sambil memberi penalti agar kurva tidak terlalu berlekuk. Tujuannya adalah mengikuti pola utama data tanpa terlalu mengikuti variasi acak pada setiap titik.

Metode ini mempertimbangkan dua hal sekaligus: seberapa dekat kurva dengan data dan seberapa besar kelengkungannya. Keseimbangan keduanya diatur oleh parameter penghalusan, yang ditulis sebagai \(\lambda\). Penalti yang lebih kuat menghasilkan kurva yang lebih halus; penalti yang lebih lemah memungkinkan kurva lebih lentur.

Sebagai contoh, pengukuran pertumbuhan tanaman dapat memperlihatkan naik-turun kecil akibat variasi pengukuran. Smoothing spline membantu menggambarkan kecenderungan pertumbuhannya melalui kurva yang mulus, tanpa mengharuskan kurva melewati setiap pengamatan.

Pada regression spline, perhatian kita terutama pada pemilihan jumlah dan lokasi knot. Pada smoothing spline, perhatian utamanya adalah tingkat penghalusan. Smoothing spline tetap menggunakan knot, tetapi kita tidak perlu memilih titik perubahan satu per satu seperti pada contoh regression spline sebelumnya. Fungsi R yang digunakan dapat memilih tingkat penghalusan dari data.

8.2 Mengatur kehalusan kurva

Pada regression spline, kita menentukan jumlah dan lokasi knot. Pada smoothing spline, perhatian utama beralih ke seberapa halus kurva yang dibentuk. Kurva dipilih dengan meminimumkan

\[ \sum_{i=1}^{n}\{y_i-g(x_i)\}^2+ \lambda\int\{g''(t)\}^2\,dt. \]

Bagian pertama mengukur selisih antara data dan kurva. Bagian kedua memberi penalti pada kelengkungan kurva. Parameter \(\lambda\) mengatur keseimbangan keduanya.

• Ketika \(\lambda\) kecil, penalti lebih lemah sehingga kurva lebih leluasa mengikuti data.
• Ketika \(\lambda\) besar, kurva menjadi lebih halus dan cenderung mendekati garis lurus.

Kurva yang sangat lentur dapat mengikuti variasi acak pada data latih. Sebaliknya, kurva yang terlalu kaku dapat melewatkan pola hubungan yang perlu dijelaskan.

8.3 Menduga model di R

ss1 <- smooth.spline(dat$data.x, dat$y, all.knots = TRUE)
ss2 <- smooth.spline(dat$data.x, dat$y)

ss1

## Call:
## smooth.spline(x = dat$data.x, y = dat$y, all.knots = TRUE)
## 
## Smoothing Parameter  spar= 1.5  lambda= 0.0002746 (26 iterations)
## Equivalent Degrees of Freedom (Df): 14.85
## Penalized Criterion (RSS): 24848
## GCV: 25.6

ss2

## Call:
## smooth.spline(x = dat$data.x, y = dat$y)
## 
## Smoothing Parameter  spar= 1.108  lambda= 0.01534 (11 iterations)
## Equivalent Degrees of Freedom (Df): 5.86
## Penalized Criterion (RSS): 25104
## GCV: 25.4

Pada ss1, setiap nilai unik \(x\) dipakai sebagai knot. Pada ss2, R menggunakan pengaturan bawaan yang dapat mengurangi jumlah knot ketika data cukup besar. Secara bawaan, tingkat penghalusan dipilih menggunakan generalized cross-validation atau GCV.

Perhatikan bagian Equivalent Degrees of Freedom pada keluaran. Derajat bebas efektif menggambarkan kelenturan kurva: nilai yang lebih besar menunjukkan kurva yang lebih lentur. Garis lurus dengan intersep memiliki dua derajat bebas.

plot_sebaran(dat$data.x, dat$y, "Hasil smoothing spline")
lines(ss1, col = pal["smoothing"], lwd = 3)
lines(ss2, col = pal["spline"], lwd = 3, lty = 2)
curve(
  5 + 2 * x + 3 * x^2,
  from = min(dat$data.x), to = max(dat$data.x),
  add = TRUE, col = pal["asli"], lwd = 2, lty = 3
)
legend(
  "topleft",
  c("Seluruh knot", "Pengaturan bawaan", "Kurva pembangkit data"),
  col = pal[c("smoothing", "spline", "asli")],
  lwd = c(3, 3, 2), lty = c(1, 2, 3), bty = "n", cex = 0.85
)

Garis hitam menunjukkan kurva yang kita gunakan untuk membangkitkan data. Garis ini dapat ditampilkan karena kita menggunakan simulasi. Pada data nyata, bentuk hubungan sebenarnya tidak diketahui.

Perhatikan juga ujung kiri dan kanan grafik. Jika pengamatan pada bagian tersebut sedikit, bentuk kurva lebih sulit diduga dengan stabil. Karena itu, pembacaan kurva perlu mempertimbangkan sebaran titik-titiknya.

8.4 Melihat pengaruh derajat bebas

Kita dapat mengatur derajat bebas untuk melihat perbedaan kelenturan kurva.

ss_kaku <- smooth.spline(dat$data.x, dat$y, df = 3)
ss_lentur <- smooth.spline(dat$data.x, dat$y, df = 40)

plot_sebaran(dat$data.x, dat$y, "Perbandingan kelenturan smoothing spline")
lines(ss_kaku, col = pal["linear"], lwd = 3)
lines(ss2, col = pal["smoothing"], lwd = 3)
lines(ss_lentur, col = pal["tangga"], lwd = 2)

legend(
  "topleft",
  c("df = 3", paste0("df otomatis = ", round(ss2$df, 1)), "df = 40"),
  col = pal[c("linear", "smoothing", "tangga")],
  lwd = c(3, 3, 2), bty = "n"
)

Amati bagian kurva yang berubah ketika derajat bebas dinaikkan. Apakah lekukan tambahan mengikuti pola utama atau hanya mengikuti variasi titik di sekitarnya?

Kurva dengan derajat bebas besar lebih berisiko mengalami overfitting, sedangkan kurva yang terlalu kaku berisiko mengalami underfitting. Namun, nilai derajat bebas saja belum cukup untuk memutuskan bahwa suatu model mengalami masalah tersebut. Kita perlu melihat kemampuan prediksinya.

9 Membandingkan model dengan cross-validation

9.1 Mengapa data perlu dibagi?

Model yang lebih rumit biasanya dapat lebih dekat dengan data latih. Namun, yang ingin kita ketahui adalah kemampuannya memprediksi respons pada pengamatan yang tidak dipakai saat pendugaan.

Pada 5-fold cross-validation, data dibagi menjadi lima bagian. Satu bagian digunakan sebagai data uji dan empat bagian lainnya sebagai data latih. Proses diulang sampai setiap bagian pernah menjadi data uji.

Galat prediksi diukur menggunakan mean squared error:

\[ \operatorname{MSE}=\frac{1}{n_{\mathrm{uji}}} \sum_{i\in\mathrm{uji}}(y_i-\widehat{y}_i)^2. \]

Semakin kecil MSE, semakin kecil rata-rata kuadrat galat prediksinya. Akar MSE disebut RMSE dan memiliki satuan yang sama dengan respons.

9.2 Menyiapkan pembagian yang sama untuk semua model

Semua model di bagian ini menggunakan dat, yaitu data dengan simpangan baku galat 5. Kita juga menggunakan pembagian fold yang sama agar perbandingannya konsisten.

set.seed(456)
fold_id <- sample(rep(1:5, length.out = nrow(dat)))
table(fold_id)

## fold_id
##   1   2   3   4   5 
## 200 200 200 200 200

Kode tersebut membuat lima kelompok dengan ukuran seimbang, lalu mengacak pengamatan yang masuk ke setiap kelompok. Pembagian acak sesuai untuk contoh simulasi ini yang pengamatannya dibangkitkan secara independen. Pada data berurutan menurut waktu atau memiliki pengelompokan, pembagian validasi perlu disesuaikan.

9.3 Membuat fungsi evaluasi

Fungsi berikut menduga model lm() pada data latih, membuat prediksi pada data uji, lalu menghitung MSE untuk setiap fold.

cv_mse <- function(formula, data, fold = fold_id) {
  hasil <- numeric(length(unique(fold)))

  for (i in sort(unique(fold))) {
    latih <- data[fold != i, , drop = FALSE]
    uji <- data[fold == i, , drop = FALSE]

    model <- lm(formula, data = latih)
    prediksi <- predict(model, newdata = uji)
    hasil[i] <- mean((uji$y - prediksi)^2)
  }

  hasil
}

Untuk smoothing spline, cara mengambil nilai prediksinya sedikit berbeda, sehingga kita buat fungsi tersendiri.

cv_ss <- function(data, fold = fold_id) {
  hasil <- numeric(length(unique(fold)))

  for (i in sort(unique(fold))) {
    latih <- data[fold != i, , drop = FALSE]
    uji <- data[fold == i, , drop = FALSE]

    model <- smooth.spline(latih$data.x, latih$y)
    prediksi <- predict(model, x = uji$data.x)$y
    hasil[i] <- mean((uji$y - prediksi)^2)
  }

  hasil
}

Tingkat penghalusan smoothing spline dipilih menggunakan data latih pada setiap pengulangan. Data uji tidak dipakai untuk memilih tingkat penghalusan.

9.4 Membandingkan satu dan dua knot

f_spline1 <- y ~ data.x + I(data.x^2) + I(data.x^3) +
  I(pmax(data.x - 1, 0)^3)

f_spline2 <- y ~ data.x + I(data.x^2) + I(data.x^3) +
  I(pmax(data.x, 0)^3) + I(pmax(data.x - 2, 0)^3)

cv1 <- cv_mse(f_spline1, dat)
cv2 <- cv_mse(f_spline2, dat)

knitr::kable(
  data.frame(
    Fold = c(as.character(1:5), "Rata-rata"),
    Satu_knot = c(cv1, mean(cv1)),
    Dua_knot = c(cv2, mean(cv2))
  ),
  digits = 3, row.names = FALSE
)

Fold	Satu_knot	Dua_knot
1	27.13	27.17
2	24.63	24.59
3	22.48	22.48
4	26.90	26.97
5	26.32	26.55
Rata-rata	25.49	25.55

Bandingkan baris rata-rata dan perhatikan pula variasi hasil antar-fold. Jika selisih rata-ratanya kecil, belum ada alasan kuat untuk menganggap penambahan knot memberikan perbaikan prediksi yang berarti. Model yang lebih sederhana dapat menjadi pilihan yang masuk akal.

9.5 Membandingkan seluruh metode

cv_linear <- cv_mse(y ~ data.x, dat)
cv_poli <- cv_mse(y ~ data.x + I(data.x^2), dat)
cv_tangga <- cv_mse(
  y ~ cut(data.x, breaks = c(-Inf, 0, 2, Inf)),
  dat
)
cv_smooth <- cv_ss(dat)

mse_fold <- cbind(
  Linear = cv_linear,
  Polinomial = cv_poli,
  Tangga = cv_tangga,
  Spline_1_knot = cv1,
  Spline_2_knot = cv2,
  Smoothing_spline = cv_smooth
)

perbandingan <- data.frame(
  Model = c(
    "Regresi linear", "Polinomial derajat dua", "Fungsi tangga",
    "Spline kubik: satu knot", "Spline kubik: dua knot",
    "Smoothing spline"
  ),
  MSE_CV = colMeans(mse_fold),
  row.names = NULL
)
perbandingan$RMSE_CV <- sqrt(perbandingan$MSE_CV)
perbandingan <- perbandingan[order(perbandingan$MSE_CV), ]

knitr::kable(perbandingan, digits = 3, row.names = FALSE)

Model	MSE_CV	RMSE_CV
Smoothing spline	25.47	5.047
Spline kubik: satu knot	25.49	5.049
Polinomial derajat dua	25.53	5.053
Spline kubik: dua knot	25.55	5.055
Regresi linear	41.77	6.463
Fungsi tangga	53.83	7.337

Tabel diurutkan dari rata-rata MSE terkecil. Karena ukuran kelima fold sama, rata-rata MSE antar-fold juga sama dengan rata-rata kuadrat galat seluruh prediksi uji.

Hubungkan hasilnya dengan cara data dibentuk. Hubungan sebenarnya kuadratik, sehingga polinomial derajat dua memang memiliki bentuk yang sesuai. Spline juga dapat mendekati hubungan tersebut, sedangkan garis lurus dan tiga interval fungsi tangga memiliki keterbatasan bentuk.

Ragam galat simulasi adalah \(5^2=25\). Jika pola rata-rata berhasil ditangkap, MSE prediksi diharapkan mendekati angka ini. MSE empiris tidak harus tepat 25, bahkan bisa sedikit lebih kecil karena variasi sampel. Angka 25 merupakan ragam galat pada proses pembangkitan data, bukan batas mutlak hasil evaluasi pada setiap sampel.

Selisih MSE yang sangat kecil dalam satu pembagian cross-validation belum cukup untuk menyatakan satu metode selalu lebih baik. Hasil simulasi ini juga tidak menentukan metode terbaik untuk semua data.

9.6 Melihat seluruh kurva pada data yang sama

Model-model berikut diduga ulang pada seluruh dat untuk menggambar kurva akhirnya. Nilai evaluasi pada tabel sebelumnya tetap berasal dari prediksi data uji.

m_lin <- lm(y ~ data.x, data = dat)
m_pol <- lm(y ~ data.x + I(data.x^2), data = dat)
m_step <- lm(y ~ cut(data.x, breaks = c(-Inf, 0, 2, Inf)), data = dat)
m_spl <- lm(f_spline1, data = dat)
m_ss <- smooth.spline(dat$data.x, dat$y)

xg <- seq(min(dat$data.x), max(dat$data.x), length.out = 400)
nd <- data.frame(data.x = xg)

plot_sebaran(dat$data.x, dat$y, "Perbandingan kurva pada data yang sama")
lines(xg, predict(m_lin, nd), col = pal["linear"], lwd = 2)
lines(xg, predict(m_pol, nd), col = pal["polinomial"], lwd = 3)
lines(xg, predict(m_spl, nd), col = pal["spline"], lwd = 2, lty = 2)
lines(xg, predict(m_ss, x = xg)$y,
      col = pal["smoothing"], lwd = 2, lty = 3)

batas <- c(min(dat$data.x), 0, 2, max(dat$data.x))
tengah <- (batas[1:3] + batas[2:4]) / 2
pred_tangga <- predict(m_step, data.frame(data.x = tengah))
segments(batas[1:3], pred_tangga, batas[2:4], pred_tangga,
         col = pal["tangga"], lwd = 2)

lines(xg, 5 + 2 * xg + 3 * xg^2,
      col = pal["asli"], lwd = 2, lty = 4)

legend(
  "topleft", ncol = 2, bty = "n", cex = 0.8,
  legend = c("Linear", "Polinomial", "Fungsi tangga",
             "Spline kubik", "Smoothing spline", "Kurva pembangkit"),
  col = pal[c("linear", "polinomial", "tangga",
              "spline", "smoothing", "asli")],
  lty = c(1, 1, 1, 2, 3, 4), lwd = 2
)

Beberapa kurva mungkin hampir berhimpit. Ini menunjukkan bahwa metode dengan cara pembentukan berbeda dapat menghasilkan prediksi yang mirip pada data tertentu.

10 Hal yang perlu diingat

Metode	Kelebihan	Hal yang perlu diperhatikan
Polinomial	Mudah dibentuk untuk lengkungan sederhana.	Derajat tinggi dapat menghasilkan kurva tidak stabil, terutama di ujung rentang data.
Fungsi tangga	Mudah dibaca sebagai perbedaan rata-rata antarinterval.	Prediksi konstan dalam interval dan bergantung pada batas kelompok.
Regression spline	Dapat mengikuti perubahan pola secara mulus.	Jumlah dan lokasi knot memengaruhi kelenturan; spline kubik biasa tetap dapat tidak stabil di ujung.
Smoothing spline	Tingkat penghalusan dapat dipilih dari data.	Kurva yang lebih lentur belum tentu lebih baik untuk prediksi.

Saat menganalisis data, mulailah dengan diagram pencar. Setelah model diduga, periksa kurva dan residualnya. Gunakan evaluasi prediksi untuk menilai apakah tambahan kerumitan model memang membantu.

Pada praktikum ini, nilai \(x\) untuk menggambar prediksi dibatasi pada rentang data. Prediksi jauh di luar rentang pengamatan perlu dipertimbangkan dengan hati-hati karena bentuk hubungan di sana tidak didukung oleh data yang tersedia.

11 Latihan mandiri

11.1 Latihan 1: Data mobil

Gunakan data mtcars. Jadikan hp sebagai peubah penjelas dan mpg sebagai respons.

dat_mt <- data.frame(
  data.x = mtcars$hp,
  y = mtcars$mpg
)

1. Buat diagram pencar dan jelaskan pola hubungannya.
2. Duga regresi linear dan polinomial derajat dua.
3. Duga spline kubik dengan satu knot di median hp. Tetapkan nilai knot tersebut sebelum membandingkan model, dan gunakan nilai yang sama pada semua fold latihan.
4. Gambarkan ketiga kurva dalam satu grafik.
5. Bandingkan MSE menggunakan 5-fold cross-validation dengan pembagian fold yang sama.
6. Jelaskan apakah model dengan MSE terkecil juga terlihat sesuai pada grafik.

Gunakan kode berikut untuk menyiapkan pembagian fold. Argumen fold harus diberikan secara eksplisit karena data latihan berbeda dari data simulasi.

set.seed(456)
fold_mt <- sample(rep(1:5, length.out = nrow(dat_mt)))

cv_mt_linear <- cv_mse(
  y ~ data.x,
  data = dat_mt,
  fold = fold_mt
)
mean(cv_mt_linear)

Jumlah pengamatan hanya 32. Karena itu, hasil evaluasi dapat cukup berubah jika pembagian fold diganti. Pertimbangkan keterbatasan ini ketika menarik kesimpulan.

11.2 Latihan 2: Mengubah knot

Gunakan kembali data simulasi dat.

1. Tentukan tiga knot berdasarkan kuartil \(x\), yaitu kuantil 0,25; 0,50; dan 0,75.
2. Tambahkan satu basis pangkat terpotong untuk setiap knot.
3. Hitung MSE menggunakan pembagian fold_id yang sudah tersedia.
4. Bandingkan dengan spline satu dan dua knot. Apakah tambahan knot memperbaiki prediksi?