Praktikum 10: Supervised dan Unsupervised Learning

1. Praktik dengan data padi

Dari data yang sama, kita bisa membuat prediksi atau mencari kelompok. Hari ini kita mencoba keduanya lewat studi kasus pertanian.

Pertanyaan	Pendekatan	Hasil yang dicari
Berapa hasil panen pada kondisi tertentu?	Supervised learning: regresi	Prediksi angka
Kondisi lahan mana yang mirip?	Unsupervised learning: clustering	Kelompok pengamatan

Pada akhir praktikum, kita bisa membaca pola data, menambah fitur yang sesuai, menilai prediksi pada data uji, dan menjelaskan hasil K-Means. Waktu praktik sekitar dua jam.

Sebelum mulai: pada pertanyaan pertama, kolom mana yang menjadi target? Apakah kolom itu juga harus dipakai pada pertanyaan kedua?

Mulai di R

Ganti nim_seed sesuai arahan kelas. Jalankan kode secara berurutan; tombol Code membuka sintaks di setiap bagian.

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

nim_seed <- 123
set.seed(nim_seed)
n <- 150

padi <- data.frame(
  curah_hujan = rnorm(n, mean = 180, sd = 40),
  suhu = rnorm(n, mean = 27, sd = 2),
  pupuk = rnorm(n, mean = 110, sd = 20)
)

padi$hasil_panen <- -0.0008 * (padi$curah_hujan - 180)^2 +
  0.03 * padi$pupuk + rnorm(n, mean = 0, sd = 1)

knitr::kable(head(padi, 5), digits = 2,
  col.names = c("Curah hujan", "Suhu", "Pupuk", "Skor hasil panen"),
  row.names = FALSE)

Curah hujan	Suhu	Pupuk	Skor hasil panen
157.6	28.58	95.70	3.90
170.8	28.54	94.95	3.83
242.3	27.66	91.23	0.06
182.8	24.98	88.95	3.38
185.2	26.76	101.26	3.93

Ini adalah data simulasi, mengikuti mekanisme pada modul praktikum. hasil_panen dibaca sebagai skor simulasi, bukan ton per hektare. Nilai negatif tidak berarti produksi fisik yang negatif.

Variabel	Peran pada regresi	Peran pada clustering hari ini
`curah_hujan`	Prediktor	Dasar pengelompokan
`suhu`	Prediktor	Tidak digunakan
`pupuk`	Prediktor	Dasar pengelompokan
`hasil_panen`	Target yang diprediksi	Tidak digunakan

2. Data latih dan data uji

Model belajar dari data latih. Data uji dipakai untuk memeriksa prediksi pada pengamatan yang belum ikut melatih model.

set.seed(nim_seed + 1)
id_latih <- sample(seq_len(nrow(padi)), size = floor(0.8 * nrow(padi)))

latih <- padi[id_latih, ]
uji <- padi[-id_latih, ]

knitr::kable(data.frame(
  Bagian = c("Data latih", "Data uji"),
  Jumlah = c(nrow(latih), nrow(uji)),
  Kegunaan = c("EDA dan pelatihan model", "Evaluasi prediksi")
), row.names = FALSE)

Bagian	Jumlah	Kegunaan
Data latih	120	EDA dan pelatihan model
Data uji	30	Evaluasi prediksi

knitr::kable(data.frame(
  Variabel = names(latih),
  Nilai_hilang = colSums(is.na(latih))
), col.names = c("Variabel", "Nilai hilang"), row.names = FALSE)

Variabel	Nilai hilang
curah_hujan	0
suhu	0
pupuk	0
hasil_panen	0

EDA untuk membentuk model berikut memakai data latih. Dengan begitu, kita tidak mengintip pola target pada data uji untuk merancang fitur.

Diskusi singkat: kalau model dilatih dan dinilai pada pengamatan yang sama, kemampuan apa yang sebenarnya sedang kita ukur?

3. EDA sebelum pemodelan

Distribusi dan pengamatan ekstrem

Mulai dari dua hal: rentang curah hujan dan sebaran skor hasil panen.

par(mfrow = c(1, 2))
tema_plot()
hist(latih$curah_hujan, breaks = 12,
     col = "#AAD1DF", border = "white",
     main = "Sebaran curah hujan", xlab = "Curah hujan", ylab = "Frekuensi")

tema_plot()
boxplot(latih$hasil_panen, horizontal = TRUE,
        col = "#F0C59E", border = jingga,
        main = "Sebaran skor hasil panen", xlab = "Skor hasil panen", yaxt = "n")
stripchart(latih$hasil_panen, method = "jitter", add = TRUE,
           vertical = FALSE, pch = 16, cex = 0.7,
           col = adjustcolor(biru, alpha.f = 0.5))

par(mfrow = c(1, 1))

Apakah ada pengamatan yang jauh dari sebagian besar data?
Apakah nilai ekstrem itu kesalahan pencatatan atau bagian dari pola yang perlu dipelajari?

Titik di luar whisker boxplot tidak otomatis salah. Periksa konteksnya sebelum memutuskan untuk menghapusnya. Prediktor juga tidak harus berdistribusi normal agar dapat masuk ke regresi linear.

Hubungan prediktor dengan hasil panen

par(mfrow = c(1, 3))

tema_plot()
plot(latih$curah_hujan, latih$hasil_panen,
     pch = 19, cex = 0.85, col = adjustcolor(biru, 0.6),
     xlab = "Curah hujan", ylab = "Skor hasil panen", main = "Curah hujan")
lines(lowess(latih$curah_hujan, latih$hasil_panen), col = jingga, lwd = 3)

tema_plot()
plot(latih$suhu, latih$hasil_panen,
     pch = 19, cex = 0.85, col = adjustcolor(biru, 0.6),
     xlab = "Suhu", ylab = "Skor hasil panen", main = "Suhu")
lines(lowess(latih$suhu, latih$hasil_panen), col = jingga, lwd = 3)

tema_plot()
plot(latih$pupuk, latih$hasil_panen,
     pch = 19, cex = 0.85, col = adjustcolor(biru, 0.6),
     xlab = "Pupuk", ylab = "Skor hasil panen", main = "Pupuk")
lines(lowess(latih$pupuk, latih$hasil_panen), col = jingga, lwd = 3)

par(mfrow = c(1, 1))

Titik adalah pengamatan pada data latih. Garis jingga merangkum kecenderungan lokal dengan lowess(); garis ini membantu EDA, belum menjadi model prediksi yang akan kita evaluasi.

korelasi <- cor(latih[c("curah_hujan", "suhu", "pupuk")], latih$hasil_panen)
knitr::kable(data.frame(
  Prediktor = rownames(korelasi),
  Korelasi_linear = as.numeric(korelasi)
), digits = 3, col.names = c("Prediktor", "Korelasi linear dengan target"),
row.names = FALSE)

Prediktor	Korelasi linear dengan target
curah_hujan	-0.213
suhu	0.094
pupuk	0.377

Amati grafik curah hujan: apakah satu garis lurus sudah cukup? Bandingkan bentuk grafiknya dengan nilai korelasi linear pada tabel.

Korelasi linear yang kecil dapat muncul meskipun hubungan melengkung terlihat jelas. Pola pada simulasi ini memberi alasan untuk mencoba fitur kuadratik. Pola tersebut tidak boleh langsung ditafsirkan sebagai bukti pengaruh sebab akibat pada lahan nyata.

4. Fitur untuk pola melengkung

Kita tambahkan fitur yang mengukur kuadrat selisih curah hujan dari 180. Nilai 180 sudah ditetapkan dalam skenario simulasi, bukan dicari dari data uji.

latih$hujan_kuadrat <- (latih$curah_hujan - 180)^2
uji$hujan_kuadrat <- (uji$curah_hujan - 180)^2

contoh_hujan <- c(130, 155, 180, 205, 230)
knitr::kable(data.frame(
  Curah_hujan = contoh_hujan,
  Selisih_dari_180 = contoh_hujan - 180,
  Kuadrat_selisih = (contoh_hujan - 180)^2
), col.names = c("Curah hujan", "Selisih dari 180", "Kuadrat selisih"),
row.names = FALSE)

Curah hujan	Selisih dari 180	Kuadrat selisih
130	-50	2500
155	-25	625
180	0	0
205	25	625
230	50	2500

Curah hujan yang berjarak sama dari 180 menghasilkan nilai kuadrat yang sama. Fitur ini memberi ruang bagi model untuk membentuk pola melengkung.

ŷ = β₀ + β₁ hujan + β₂ (hujan − 180)² + β₃ suhu + β₄ pupuk

Hubungannya boleh melengkung, tetapi model ini tetap linear terhadap parameter β. Kita masih bisa menduganya menggunakan lm().

5. Regresi untuk prediksi

Dua bentuk model

Kita membandingkan spesifikasi dengan prediktor linear dan spesifikasi dengan tambahan fitur kuadratik. Keduanya memakai data latih yang sama.

model_linear <- lm(
  hasil_panen ~ curah_hujan + suhu + pupuk,
  data = latih
)

model_kuadratik <- lm(
  hasil_panen ~ curah_hujan + hujan_kuadrat + suhu + pupuk,
  data = latih
)

knitr::kable(data.frame(
  Komponen = names(coef(model_kuadratik)),
  Koefisien = unname(coef(model_kuadratik))
), digits = 5, row.names = FALSE)

Komponen	Koefisien
(Intercept)	0.28615
curah_hujan	-0.00069
hujan_kuadrat	-0.00086
suhu	-0.00275
pupuk	0.03062

Perhatikan tanda koefisien hujan_kuadrat. Nilai negatif memungkinkan pola melengkung ke bawah. Besar koefisien mentah antarvariabel tidak langsung menunjukkan variabel yang paling penting karena skala dan bentuk fiturnya berbeda.

Bentuk prediksi yang dipelajari

grid_hujan <- data.frame(
  curah_hujan = seq(min(latih$curah_hujan), max(latih$curah_hujan), length.out = 180),
  suhu = mean(latih$suhu),
  pupuk = mean(latih$pupuk)
)
grid_hujan$hujan_kuadrat <- (grid_hujan$curah_hujan - 180)^2

tema_plot()
plot(latih$curah_hujan, latih$hasil_panen,
     pch = 19, col = adjustcolor("#5D7E8D", 0.45),
     xlab = "Curah hujan", ylab = "Skor hasil panen",
     main = "Garis lurus dan pola kuadratik")
lines(grid_hujan$curah_hujan, predict(model_linear, grid_hujan),
      col = jingga, lwd = 3, lty = 2)
lines(grid_hujan$curah_hujan, predict(model_kuadratik, grid_hujan),
      col = biru, lwd = 3)
legend("bottom", horiz = TRUE, bty = "n",
       legend = c("Linear", "Dengan fitur kuadratik"),
       col = c(jingga, biru), lwd = 3, lty = c(2, 1), cex = 0.9)

Kurva prediksi menahan suhu dan pupuk pada rata-rata data latih. Titik tetap memakai nilai suhu dan pupuk masing-masing, sehingga tidak harus tepat berada pada kurva.

Prediksi pada data uji

pred_linear <- predict(model_linear, newdata = uji)
pred_kuadratik <- predict(model_kuadratik, newdata = uji)

knitr::kable(head(data.frame(
  Aktual = uji$hasil_panen,
  Prediksi_linear = pred_linear,
  Prediksi_kuadratik = pred_kuadratik
), 5), digits = 3, row.names = FALSE)

Aktual	Prediksi_linear	Prediksi_kuadratik
3.899	1.999	2.595
2.930	0.814	1.912
1.601	2.026	2.831
-0.681	1.764	-2.598
2.777	3.217	3.135

newdata = uji memastikan prediksi dihitung untuk pengamatan uji. Jika newdata dihilangkan, predict() pada model lm mengembalikan prediksi untuk data yang dipakai melatih model.

6. Evaluasi pada data uji

Error menunjukkan seberapa jauh prediksi dari nilai aktual. Kita gunakan MAE dan RMSE, keduanya dalam satuan skor hasil panen.

Ukuran	Cara membaca
MAE	Rata-rata besar selisih absolut antara aktual dan prediksi
RMSE	Akar rata-rata kuadrat error; lebih peka terhadap error yang besar

mae <- function(aktual, prediksi) mean(abs(aktual - prediksi))
rmse <- function(aktual, prediksi) sqrt(mean((aktual - prediksi)^2))

pred_baseline <- rep(mean(latih$hasil_panen), nrow(uji))

evaluasi <- data.frame(
  Model = c("Rata-rata data latih", "Linear", "Dengan fitur kuadratik"),
  MAE_latih = c(
    mae(latih$hasil_panen, rep(mean(latih$hasil_panen), nrow(latih))),
    mae(latih$hasil_panen, fitted(model_linear)),
    mae(latih$hasil_panen, fitted(model_kuadratik))
  ),
  MAE_uji = c(mae(uji$hasil_panen, pred_baseline),
              mae(uji$hasil_panen, pred_linear),
              mae(uji$hasil_panen, pred_kuadratik)),
  RMSE_uji = c(rmse(uji$hasil_panen, pred_baseline),
               rmse(uji$hasil_panen, pred_linear),
               rmse(uji$hasil_panen, pred_kuadratik))
)

knitr::kable(evaluasi, digits = 3,
  col.names = c("Model", "MAE latih", "MAE uji", "RMSE uji"), row.names = FALSE)

Model	MAE latih	MAE uji	RMSE uji
Rata-rata data latih	1.518	1.423	1.987
Linear	1.394	1.418	1.920
Dengan fitur kuadratik	0.747	0.745	0.890

MAE uji model kuadratik adalah 0.745. Artinya, pada data uji ini, rata-rata besar selisih prediksi dari skor aktual sekitar 0.745 poin.

Baseline selalu menebak rata-rata target pada data latih. Model yang lebih rumit perlu memberikan manfaat dibandingkan tebakan sederhana ini.

rentang <- range(c(uji$hasil_panen, pred_linear, pred_kuadratik))
par(mfrow = c(1, 2))

tema_plot()
plot(uji$hasil_panen, pred_linear,
     xlim = rentang, ylim = rentang, pch = 19, col = adjustcolor(jingga, 0.8),
     xlab = "Skor aktual", ylab = "Skor prediksi", main = "Model linear")
abline(a = 0, b = 1, col = "#80919B", lty = 2, lwd = 2)

tema_plot()
plot(uji$hasil_panen, pred_kuadratik,
     xlim = rentang, ylim = rentang, pch = 19, col = adjustcolor(biru, 0.8),
     xlab = "Skor aktual", ylab = "Skor prediksi", main = "Model kuadratik")
abline(a = 0, b = 1, col = "#80919B", lty = 2, lwd = 2)

par(mfrow = c(1, 1))

Semua titik pada kedua panel berasal dari data uji yang sama. Makin dekat titik ke garis putus-putus, makin dekat prediksi ke nilai aktual.

Diskusi: apakah error harus menjadi nol agar model berguna? Apa yang terjadi jika error data latih sangat kecil, tetapi error data uji besar?

Perbandingan ini dipakai untuk membaca kinerja dua spesifikasi yang telah ditentukan. Jika ingin mencari model atau parameter terbaik dari banyak percobaan, lakukan validasi atau cross-validation pada data latih. Simpan data uji untuk penilaian akhir.

7. K-Means untuk pengelompokan

Sekarang pertanyaannya berubah: pengamatan mana yang memiliki kombinasi curah hujan dan pupuk yang mirip? Kita memakai seluruh pengamatan untuk eksplorasi clustering, tanpa memasukkan hasil_panen.

Samakan skala sebelum menghitung jarak

X_cluster <- padi[c("curah_hujan", "pupuk")]
X_standar <- scale(X_cluster)

knitr::kable(data.frame(
  Variabel = colnames(X_cluster),
  SD_awal = sapply(X_cluster, sd),
  Rata_rata_setelah_scale = colMeans(X_standar),
  SD_setelah_scale = apply(X_standar, 2, sd)
), digits = 3,
col.names = c("Variabel", "SD awal", "Rata-rata setelah scale", "SD setelah scale"),
row.names = FALSE)

Variabel	SD awal	Rata-rata setelah scale	SD setelah scale
curah_hujan	37.99	0	1
pupuk	20.63	0	1

Nilai standar = (nilai awal − rata-rata) / simpangan baku

scale() membuat kedua fitur memiliki rata-rata sekitar 0 dan simpangan baku 1. Ini membantu agar perbedaan skala tidak membuat satu variabel mendominasi jarak.

Contoh partisi menjadi tiga cluster

set.seed(nim_seed)
km3 <- kmeans(X_standar, centers = 3, nstart = 25, iter.max = 100)

centers = 3: banyak cluster untuk contoh ini.
nstart = 25: mencoba 25 inisialisasi dan mengambil hasil dengan SSE terkecil di antara percobaan tersebut.
iter.max = 100: batas iterasi untuk satu percobaan.

tema_plot()
plot(padi$curah_hujan, padi$pupuk,
     pch = 19, col = adjustcolor(warna_cluster[km3$cluster], 0.8),
     xlab = "Curah hujan", ylab = "Pupuk",
     main = "K-Means dengan K = 3")
pusat3 <- aggregate(X_cluster, by = list(Cluster = km3$cluster), FUN = mean)
points(pusat3$curah_hujan, pusat3$pupuk, pch = 4, cex = 2.1, lwd = 3)
legend("topright", legend = paste("Cluster", 1:3),
       col = warna_cluster[1:3], pch = 19, bty = "n", cex = 0.85)

K-Means dihitung pada data standar. Grafik ditampilkan kembali pada skala awal agar mudah dibaca. Tanda silang menunjukkan rata-rata setiap cluster pada skala awal.

profil <- pusat3
profil$Jumlah <- as.integer(table(km3$cluster))
profil <- profil[c("Cluster", "Jumlah", "curah_hujan", "pupuk")]
knitr::kable(profil, digits = 2,
  col.names = c("Cluster", "Jumlah pengamatan", "Rata-rata curah hujan", "Rata-rata pupuk"),
  row.names = FALSE)

Cluster	Jumlah pengamatan	Rata-rata curah hujan	Rata-rata pupuk
1	60	155.4	125.5
2	48	169.0	88.9
3	42	224.3	108.7

Gunakan tabel ini untuk memberi deskripsi, misalnya kelompok dengan rata-rata curah hujan lebih tinggi. Nomor cluster tidak menunjukkan urutan kualitas atau produktivitas lahan.

8. Jumlah cluster dan kualitasnya

Elbow dan silhouette

SSE mengukur jumlah kuadrat jarak ke pusat cluster. Penambahan cluster biasanya menurunkan SSE, sehingga SSE terkecil saja belum cukup untuk menentukan K.

K <- 1:6
SSE <- numeric(length(K))
siluet <- rep(NA_real_, length(K))
model_k <- vector("list", length(K))

for (k in K) {
  set.seed(nim_seed)
  model_k[[k]] <- kmeans(X_standar, centers = k, nstart = 25, iter.max = 100)
  SSE[k] <- model_k[[k]]$tot.withinss
  if (k >= 2) {
    siluet[k] <- mean(cluster::silhouette(
      model_k[[k]]$cluster, dist(X_standar)
    )[, "sil_width"])
  }
}

knitr::kable(data.frame(K = K, SSE = SSE, Silhouette = siluet),
             digits = 3, row.names = FALSE)

K	SSE	Silhouette
1	298.00	NA
2	185.90	0.334
3	127.80	0.352
4	100.84	0.330
5	79.44	0.341
6	64.39	0.348

par(mfrow = c(1, 2))
tema_plot()
plot(K, SSE, type = "b", pch = 19, lwd = 2.5, col = biru,
     xlab = "Jumlah cluster (K)", ylab = "SSE", main = "Elbow method", xaxt = "n")
axis(1, at = K)

tema_plot()
plot(K[-1], siluet[-1], type = "b", pch = 19, lwd = 2.5, col = hijau,
     xlab = "Jumlah cluster (K)", ylab = "Rata-rata silhouette",
     main = "Silhouette", xaxt = "n", ylim = c(0, 1))
axis(1, at = K[-1])

par(mfrow = c(1, 1))

Yang diperhatikan	Cara membaca
Elbow	Cari perubahan saat tambahan K mulai memberi penurunan SSE yang lebih kecil. Jika tekukannya tidak jelas, jangan dipaksakan.
Silhouette mendekati 1	Pengamatan relatif dekat dengan kelompoknya dan jauh dari kelompok lain.
Silhouette sekitar 0	Pengamatan berada dekat batas antarkelompok.
Silhouette negatif	Pengamatan dapat lebih dekat dengan kelompok lain.

Silhouette tidak didefinisikan untuk satu cluster; karena itu nilai K = 1 pada tabel ditulis NA.

Coba ganti K

Pilih jumlah cluster, lalu perhatikan titik mana yang berpindah kelompok. Grafik ini memakai hasil perhitungan R yang sama dengan tabel di atas.

Jumlah cluster

K = 3; SSE = 127.799; rata-rata silhouette = 0.352.

K = 3 di bagian sebelumnya adalah contoh untuk praktik, bukan jumlah kelompok alami yang sudah diketahui. Data ini dibangkitkan dari distribusi kontinu, tanpa membuat tiga populasi terpisah. Hasil partisi perlu dibaca bersama kurva evaluasi dan tujuan analisis.

Diskusi: kalau K dinaikkan terus, apakah hasilnya pasti lebih berguna? Apakah mengganti nama Cluster 1 menjadi Cluster 2 mengubah kualitas pengelompokannya?

9. Latihan penutup

A. Regresi dan fitur

Gunakan pembagian data yang sama. Coba model kuadratik tanpa suhu.

model_tanpa_suhu <- lm(
  hasil_panen ~ curah_hujan + hujan_kuadrat + pupuk,
  data = latih
)
pred_tanpa_suhu <- predict(model_tanpa_suhu, newdata = uji)
mae(uji$hasil_panen, pred_tanpa_suhu)

Bandingkan MAE-nya dengan model kuadratik awal. Jelaskan mengapa menambah satu prediktor tidak selalu memperbaiki prediksi pada data uji. Jangan menjadikan satu pembagian data sebagai bukti bahwa suatu fitur selalu tidak berguna.

B. Pengelompokan dan skala

set.seed(nim_seed)
km4 <- kmeans(X_standar, centers = 4, nstart = 25)
table(km4$cluster)

set.seed(nim_seed)
km_tanpa_scale <- kmeans(X_cluster, centers = 3, nstart = 25)

Bandingkan K = 3 dengan K = 4 melalui grafik interaktif dan tabel evaluasi. Lalu jelaskan mengapa hasil tanpa standardisasi dapat berbeda. Bandingkan SSE hanya jika memakai skala fitur yang sama.

Ringkasan hasil praktik

Bagian	Hasil yang dituliskan
Regresi	Tabel MAE model awal dan model tanpa suhu, disertai interpretasi singkat
Clustering	Perbandingan K = 3 dan K = 4 serta deskripsi profil kelompok
Refleksi	Peran EDA, alasan menambah fitur kuadratik, dan perbedaan prediksi dengan pengelompokan

Cek pemahaman sebelum selesai

Mengapa error tidak nol? Ada komponen acak dalam simulasi, parameter diduga dari sampel, dan model bisa belum menangkap seluruh pola.
Mengapa tetap perlu data uji? Kita ingin mengetahui kemampuan prediksi pada pengamatan yang tidak ikut melatih model.
Apakah hasil clustering dinilai dengan MAE hasil panen? Tidak. Target itu tidak menjadi acuan pembentukan cluster dalam praktik ini.
Apakah model lebih kompleks selalu lebih baik? Tidak. Periksa manfaatnya pada evaluasi yang sesuai dan perhatikan kompleksitasnya.

Dokumentasi fungsi R

Dokumentasi resmi: regresi linear, prediksi model, standardisasi, K-Means, dan silhouette.