Mulai di R

Buka file .Rmd ini di RStudio atau Posit Cloud. Jalankan kode secara berurutan; tombol Code membuka sintaks di setiap bagian. Gunakan Knit untuk membuat HTML.

nim_seed boleh diganti sesuai arahan kelas. Nilai awalnya disediakan agar hasil contoh dapat diulang.

Jalankan satu kali di Console jika paketnya belum tersedia:

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

Data latihan

Data simulasi. Data ini dibuat untuk pembelajaran, bukan hasil pengamatan lapangan. Hubungan di dalamnya tidak boleh dijadikan kesimpulan agronomis.

Ada 160 pengamatan dengan 10 prediktor. Target yield dinyatakan dalam ton/ha. Kode berikut membangkitkan data yang sama setiap kali dijalankan, sehingga praktikum tidak memerlukan unduhan dataset tambahan.

nim_seed <- 1101
set.seed(nim_seed)
n <- 160

air <- rnorm(n)
pupuk <- rnorm(n)
panas <- rnorm(n)

agri <- data.frame(
  rainfall      = 180 + 35 * air + rnorm(n, sd = 2),
  humidity      = 75 + 5 * air + rnorm(n, sd = 3),
  temperature   = 28 + 1.7 * panas,
  soil_moisture = 30 + 6 * air + rnorm(n, sd = 0.35),
  irrigation    = 60 + 10 * air + rnorm(n, sd = 6),
  fertilizer_A  = 120 + 20 * pupuk + rnorm(n, sd = 1),
  fertilizer_B  = 100 + 17 * pupuk + rnorm(n, sd = 1),
  sunlight      = 7 + 0.8 * panas + rnorm(n, sd = 0.4),
  soil_pH       = 6.2 + rnorm(n, sd = 0.35),
  pest_index    = runif(n, min = 0, max = 5)
)

agri$yield <- 5.5 +
  0.45 * (agri$rainfall - 180) / 35 +
  0.45 * (agri$soil_moisture - 30) / 6 +
  0.65 * (agri$fertilizer_A - 120) / 20 -
  0.45 * (agri$temperature - 28) / 1.7 -
  0.50 * (agri$pest_index - 2.5) / 1.5 +
  rnorm(n, sd = 0.85)

knitr::kable(
  head(agri[c("rainfall", "temperature", "fertilizer_A", "fertilizer_B", "yield")], 5),
  digits = 2,
  col.names = c("Curah hujan", "Suhu", "Pupuk A", "Pupuk B", "Hasil panen"),
  row.names = FALSE
)

Coba jawab: mana targetnya, mana prediktornya, dan mengapa ini termasuk regresi?

Pisahkan data latih dan data uji

Gunakan 70% data latih untuk membangun model dan memilih penalti. Simpan 30% data uji untuk evaluasi akhir. Pembagian acak sesuai untuk contoh ini karena pengamatan disimulasikan independen; data deret waktu memerlukan pembagian yang mengikuti urutan waktu.

set.seed(nim_seed + 1)
id_train <- sample(seq_len(nrow(agri)), size = floor(0.7 * nrow(agri)))
train <- agri[id_train, ]
test  <- agri[-id_train, ]

knitr::kable(data.frame(
  Bagian = c("Data latih", "Data uji"),
  Jumlah = c(nrow(train), nrow(test)),
  Kegunaan = c("EDA, pelatihan, dan cross-validation", "Evaluasi akhir")
), row.names = FALSE)

stopifnot(!anyNA(agri))

Apa yang terlihat dari dua prediktor ini?

par(mfrow = c(1, 2))
tema_plot()
plot(train$rainfall, train$yield,
     pch = 19, col = adjustcolor(warna["OLS"], 0.65),
     xlab = "Curah hujan (mm/bulan)", ylab = "Hasil panen (ton/ha)",
     main = "Hujan dan hasil panen")
tema_plot()
plot(train$pest_index, train$yield,
     pch = 19, col = adjustcolor(warna["Lasso"], 0.65),
     xlab = "Indeks hama", ylab = "Hasil panen (ton/ha)",
     main = "Hama dan hasil panen")

par(mfrow = c(1, 1))

Perhatikan arah hubungan, sebaran titik, dan kemungkinan pencilan. Scatter plot ini menunjukkan hubungan dua peubah; koefisien regresi berganda nanti memperhitungkan prediktor lain juga.

Prediktor mana yang membawa informasi serupa?

fitur <- setdiff(names(train), "yield")
label_fitur <- c("Hujan", "Kelemb. udara", "Suhu", "Kelemb. tanah",
                 "Irigasi", "Pupuk A", "Pupuk B", "Penyinaran",
                 "pH tanah", "Indeks hama")
x_train <- as.matrix(train[, fitur])
y_train <- train$yield
x_test  <- as.matrix(test[, fitur])
y_test  <- test$yield

kor <- cor(x_train)
p <- ncol(kor)
tema_plot()
par(mar = c(7, 8, 3, 2))
image(seq_len(p), seq_len(p), t(kor[p:1, ]),
      col = colorRampPalette(c("#176E91", "#FFFFFF", "#C75C45"))(101),
      zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "")
axis(1, at = seq_len(p), labels = label_fitur, las = 2, cex.axis = 0.85)
axis(2, at = seq_len(p), labels = rev(label_fitur), las = 2, cex.axis = 0.85)
for (i in seq_len(p)) for (j in seq_len(p)) {
  nilai <- kor[p + 1 - j, i]
  text(i, j, sprintf("%.3f", nilai), cex = 0.72,
       col = if (abs(nilai) > 0.65) "white" else "#243D50")
}
title("Korelasi prediktor pada data latih", line = 1)

Merah menunjukkan korelasi positif, biru menunjukkan korelasi negatif, dan warna mendekati putih menunjukkan korelasi mendekati nol. Lihat angkanya untuk membaca kekuatan hubungan.

Berhenti sebentar: temukan dua pasangan prediktor yang sangat berkorelasi. Jika keduanya masuk model, apakah mudah memisahkan kontribusinya?

Multikolinearitas terjadi ketika satu atau beberapa prediktor dapat dijelaskan kuat oleh kombinasi prediktor lain. Heatmap membantu melihat pasangan yang berkorelasi, tetapi belum menangkap seluruh kemungkinan hubungan multivariat.

Apakah koefisien besar berarti peubah paling penting?

Belum tentu. Satuan hujan, suhu, dan pH berbeda. Untuk membandingkan perubahan prediksi pada skala yang sama, kalikan setiap koefisien dengan simpangan baku prediktornya pada data latih.

sd_train <- apply(x_train, 2, sd)
beta_ols <- coef(model_ols)[fitur] * sd_train

tema_plot()
par(mar = c(5, 9, 3, 1))
barplot(rev(beta_ols), names.arg = rev(label_fitur), horiz = TRUE,
        las = 1, col = warna["OLS"], border = NA,
        xlim = range(c(0, beta_ols)) * 1.12,
        xlab = "Perubahan prediksi per kenaikan 1 SD prediktor (ton/ha)",
        main = "Koefisien OLS pada skala yang sebanding")
abline(v = 0, col = "#80919B")

Perhatikan koefisien pada pasangan yang sangat berkorelasi. Nilai ekstrem atau tanda yang sulit dijelaskan dapat menjadi alasan untuk memeriksa kestabilan, tetapi satu model saja belum membuktikan bahwa koefisiennya tidak stabil.

Ingat: multikolinearitas dapat memperbesar ketidakpastian koefisien. Prediksi gabungan masih bisa baik, terutama pada data baru dengan pola prediktor yang serupa.

Mengapa prediktor perlu distandardisasi?

Penalti bekerja pada angka koefisien. Mengubah curah hujan dari mm menjadi meter saja dapat mengubah besar koefisiennya. Standardisasi membuat penalti bekerja pada skala prediktor yang sebanding.

Pada praktik ini: berikan matriks prediktor dalam satuan asli. Gunakan standardize = TRUE agar glmnet menghitung standardisasi saat model dilatih, termasuk di dalam setiap fold CV. Koefisien yang dikembalikan tetap dalam satuan asli.

Dua pengaturan yang berbeda

Catatan saat membaca materi Python: argumen alpha pada sklearn.Ridge atau sklearn.Lasso berperan sebagai kekuatan penalti. Pada glmnet di R, kekuatan penalti diatur oleh lambda.

Lihat ketiga model dalam satu gambar

koef <- cbind(OLS = beta_ols, Ridge = beta_ridge, Lasso = beta_lasso)
rownames(koef) <- label_fitur
tema_plot()
par(mar = c(5, 9, 3, 1))
barplot(t(koef[nrow(koef):1, ]), beside = TRUE, horiz = TRUE,
        names.arg = rev(label_fitur), las = 1,
        xlim = range(c(0, koef)) * 1.12,
        col = warna, border = NA,
        xlab = "Perubahan prediksi per kenaikan 1 SD prediktor (ton/ha)",
        main = "Bagaimana koefisien berubah?",
        legend.text = colnames(koef),
        args.legend = list(x = "bottomright", bty = "n", cex = 0.85))
abline(v = 0, col = "#80919B")

Baru sekarang, buka data uji

rmse <- function(aktual, prediksi) sqrt(mean((aktual - prediksi)^2))
r2 <- function(aktual, prediksi) {
  1 - sum((aktual - prediksi)^2) / sum((aktual - mean(aktual))^2)
}

pred_train <- list(
  OLS = as.numeric(predict(model_ols, newdata = train)),
  Ridge = as.numeric(predict(cv_ridge, newx = x_train, s = "lambda.min")),
  Lasso = as.numeric(predict(cv_lasso, newx = x_train, s = "lambda.min"))
)
pred_test <- list(
  OLS = as.numeric(predict(model_ols, newdata = test)),
  Ridge = as.numeric(predict(cv_ridge, newx = x_test, s = "lambda.min")),
  Lasso = as.numeric(predict(cv_lasso, newx = x_test, s = "lambda.min"))
)

hasil <- data.frame(
  Model = names(pred_test),
  RMSE_latih = vapply(pred_train, function(z) rmse(y_train, z), numeric(1)),
  RMSE_uji = vapply(pred_test, function(z) rmse(y_test, z), numeric(1)),
  R2_uji = vapply(pred_test, function(z) r2(y_test, z), numeric(1)),
  Prediktor_aktif = c(
    sum(abs(coef(model_ols)[fitur]) > 1e-8),
    sum(abs(beta_ridge_raw) > 1e-8),
    sum(abs(beta_lasso_raw) > 1e-8)
  ), row.names = NULL
)
knitr::kable(hasil, digits = 3,
  col.names = c("Model", "RMSE latih", "RMSE uji", "R² uji", "Prediktor aktif"),
  row.names = FALSE)

• RMSE: lebih kecil lebih baik; satuannya ton/ha, sama seperti target.
• R² uji: semakin tinggi semakin baik; bisa negatif jika prediksi lebih buruk daripada pembanding berupa rata-rata nilai aktual data uji.
• Prediktor aktif: jumlah koefisien prediktor yang tidak nol secara numerik; intersep tidak dihitung.

Pada scatter plot, titik yang dekat dengan garis putus-putus menunjukkan prediksi yang dekat dengan nilai aktual. Grafik kanan menampilkan Ridge sebagai contoh, bukan sebagai pernyataan bahwa Ridge pasti paling baik.

Diskusikan: model mana yang RMSE ujinya paling kecil pada pembagian ini? Apakah hasilnya sejalan dengan jumlah prediktor aktif dan perubahan koefisiennya?

Jangan memaksakan Ridge atau Lasso menjadi pemenang. OLS juga bisa memiliki galat uji terkecil. Hasil satu pembagian data merupakan ilustrasi; selisih kecil belum cukup untuk menyatakan keunggulan yang berlaku umum. Jika data uji dipakai berulang kali untuk mengubah pengaturan, data itu sudah menjadi bagian dari proses pemilihan model.

Tambahan: benarkah koefisien lebih stabil?

Bagian ini dapat didemonstrasikan pengajar jika waktu memungkinkan. Kita mengambil ulang baris data latih sebanyak 80 kali, lalu membandingkan variasi koefisien OLS dan Ridge. Lambda Ridge ditahan pada nilai pilihan CV agar fokusnya pada perubahan sampel.

set.seed(nim_seed + 3)
B <- 80
pilih <- c("rainfall", "soil_moisture", "fertilizer_A", "fertilizer_B")
boot_ols <- boot_ridge <- matrix(NA_real_, B, length(pilih))
for (b in seq_len(B)) {
  id <- sample(seq_len(nrow(train)), replace = TRUE)
  xb <- x_train[id, , drop = FALSE]
  yb <- y_train[id]
  fit_o <- lm.fit(cbind(1, xb), yb)
  fit_r <- glmnet(xb, yb, alpha = 0, lambda = lambda_grid,
                  standardize = TRUE)
  bo <- setNames(fit_o$coefficients[-1], fitur)
  br <- as.matrix(coef(fit_r, s = cv_ridge$lambda.min))[fitur, 1]
  boot_ols[b, ] <- bo[pilih] * sd_train[pilih]
  boot_ridge[b, ] <- br[pilih] * sd_train[pilih]
}

Kotak yang lebih sempit menunjukkan variasi koefisien yang lebih kecil dalam simulasi pengambilan ulang ini. Kestabilan koefisien tidak otomatis sama dengan RMSE uji terbaik. Demonstrasi ini bersyarat pada data latih dan lambda yang sudah dipilih; ketidakpastian pemilihan lambda belum dimasukkan.

Ubah aturan pemilihan lambda Lasso

Bandingkan lambda.min dengan lambda.1se menggunakan hasil CV yang sama. Cukup periksa lambda, galat CV, dan jumlah prediktor aktif; simpan data uji untuk evaluasi akhir model yang aturannya sudah ditetapkan.

aturan <- c("lambda.min", "lambda.1se")
latihan <- do.call(rbind, lapply(aturan, function(s) {
  lam <- cv_lasso[[s]]
  idx <- which.min(abs(cv_lasso$lambda - lam))
  b <- as.matrix(coef(cv_lasso, s = s))[-1, 1]
  data.frame(Aturan = s, Lambda = lam,
             MSE_CV = cv_lasso$cvm[idx],
             Prediktor_aktif = sum(abs(b) > 1e-8))
}))
knitr::kable(latihan, digits = 4,
  col.names = c("Aturan", "Lambda", "MSE CV", "Prediktor aktif"),
  row.names = FALSE)

Kerjakan berpasangan. Siapkan jawaban singkat untuk empat pertanyaan berikut.

1. Sebutkan dua pasangan prediktor yang berkorelasi kuat dan jelaskan masalahnya bagi interpretasi koefisien OLS.
2. Bandingkan lambda.min dan lambda.1se: apa yang berubah pada penalti, galat CV, dan jumlah prediktor aktif?
3. Berdasarkan tabel evaluasi utama, model mana yang memiliki RMSE uji terkecil? Sertakan angkanya dan batasi kesimpulan pada pembagian data ini.
4. Apakah peubah yang koefisien Lasso-nya nol boleh langsung disebut tidak berpengaruh terhadap hasil panen? Jelaskan.

Luaran sederhana: file Rmd yang dapat dijalankan dan jawaban singkat atas empat pertanyaan. Tabel serta grafik utama sudah dihasilkan dari chunk praktikum.