Projek Teknik Peubah Ganda
Luthfia Aiman Hanin (G1401231064)
library(tidyverse)
library(janitor)
library(lubridate)
library(cluster)
library(factoextra)
library(reshape2)
library(scales)
library(readr)
library(rmarkdown)
library(DT)Import Data
df <- read_delim("C:/Users/lizam/Downloads/0_df_long.csv", delim=";", col_types = cols(.default = "c"))# Bersihkan nama kolom
df <- df %>% clean_names()
# Cek nama kolom
names(df)## [1] "kab_kot" "tahun" "bulan" "produk" "harga" "kategori"
## [7] "kode_bps" "kode_prov" "nama_prov"head(df, 10)## # A tibble: 10 × 9
## kab_kot tahun bulan produk harga kategori kode_bps kode_prov nama_prov
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Kab. Aceh Bar… 2022 Janu… Beras… 11429 Beras 11.05 11 Aceh
## 2 Kab. Aceh Bar… 2022 Janu… Beras… 9979 Beras 11.05 11 Aceh
## 3 Kab. Aceh Bar… 2022 Janu… Bawan… 31000 Bawang 11.05 11 Aceh
## 4 Kab. Aceh Bar… 2022 Janu… Bawan… 28636 Bawang 11.05 11 Aceh
## 5 Kab. Aceh Bar… 2022 Janu… Cabai… 19409 Cabai 11.05 11 Aceh
## 6 Kab. Aceh Bar… 2022 Janu… Cabai… 45706 Cabai 11.05 11 Aceh
## 7 Kab. Aceh Bar… 2022 Janu… Dagin… 1438… Protein 11.05 11 Aceh
## 8 Kab. Aceh Bar… 2022 Janu… Dagin… 32955 Protein 11.05 11 Aceh
## 9 Kab. Aceh Bar… 2022 Janu… Telur… 26136 Protein 11.05 11 Aceh
## 10 Kab. Aceh Bar… 2022 Janu… Gula … 14545 Gula 11.05 11 Acehn_obs <- nrow(df)
time_range <- range(paste(df$tahun, df$bulan, sep="-"))
pct_missing_before <- mean(is.na(df$harga)) * 100
cat("Jumlah observasi:", n_obs, "\n")## Jumlah observasi: 340620cat("Rentang waktu:", time_range[1], "sampai", time_range[2], "\n")## Rentang waktu: 2022-Agustus sampai 2025-Septembercat(sprintf("Persentase missing pada kolom Harga sebelum imputasi: %.2f%%\n", pct_missing_before))## Persentase missing pada kolom Harga sebelum imputasi: 18.06%Dataset berisi 340620 observasi yang mencakup rentang waktu dari Tahun 2022 Bulan Agustus sampai Tahun 2025 Bulan September. Sebelum imputasi terdapat 18.06% nilai hilang pada kolom Harga. Imputasi dilakukan secara bertingkat:
- Rata-rata pada level Produk × Kab/Kota
- Jika tidak tersedia gunakan rata-rata Produk × Provinsi
- Jika masih kosong gunakan rata-rata Produk nasional
- Fallback ke mean global.
Harga pada kabupaten/kota yang berdekatan cenderung mirip akibat biaya logistik, akses pasar, dan pola permintaan lokal sehingga memanfaatkan agregat wilayah terdekat meminimalkan bias yang ditimbulkan bila menggunakan satu nilai global saja.
# Daftar provinsi Pulau Sumatera
sumatera_prov <- c("Aceh","Sumatera Utara","Sumatera Barat","Riau","Jambi",
"Bengkulu","Sumatera Selatan","Lampung","Kepulauan Bangka Belitung","Kepulauan Riau")
unique(df$nama_prov) %>% head(40)## [1] "Aceh" "Sumatera Barat"
## [3] "Nusa Tenggara Timur" "Sumatera Utara"
## [5] "Bali" "Kalimantan Selatan"
## [7] "Jawa Barat" "Sulawesi Tengah"
## [9] "Kepulauan Bangka Belitung" "Jawa Timur"
## [11] "Jawa Tengah" "Sulawesi Selatan"
## [13] "DI Yogyakarta" "Sumatera Selatan"
## [15] "Kalimantan Tengah" "Jambi"
## [17] "Riau" "Kalimantan Barat"
## [19] "Bengkulu" "Kalimantan Timur"
## [21] "Papua Barat" "Nusa Tenggara Barat"
## [23] "Kepulauan Riau" "Gorontalo"
## [25] "Sulawesi Utara" "Sulawesi Tenggara"
## [27] NA "Maluku"
## [29] "Papua" "Maluku Utara"
## [31] "DKI Jakarta" "Lampung"
## [33] "Banten" "Sulawesi Barat"# Filter
df_ds <- df %>%
filter(produk == "Daging Sapi Murni", nama_prov %in% sumatera_prov)# ubah harga jadi numeric
df_ds <- df_ds %>%
mutate(harga = parse_number(harga))
sum(is.na(df_ds$harga)) # jumlah missing sebelum imputasi## [1] 555Imputasi Bertingkat (Wilayah Terdekat)
# Hitung rata2 pada tiap level
mean_prod_kab <- df_ds %>%
group_by(produk, kab_kot) %>%
summarise(mean_pk = mean(harga, na.rm=TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'produk'. You can override using the
## `.groups` argument.mean_prod_prov <- df_ds %>%
group_by(produk, nama_prov) %>%
summarise(mean_pp = mean(harga, na.rm=TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'produk'. You can override using the
## `.groups` argument.mean_prod_global <- df_ds %>%
group_by(produk) %>%
summarise(mean_p = mean(harga, na.rm=TRUE)) %>%
ungroup()
global_mean <- mean(df_ds$harga, na.rm=TRUE)
# Gabungkan dan imputasi
df_imp <- df_ds %>%
left_join(mean_prod_kab, by = c("produk","kab_kot")) %>%
left_join(mean_prod_prov, by = c("produk","nama_prov")) %>%
left_join(mean_prod_global, by = "produk") %>%
mutate(harga_imputed = case_when(
!is.na(harga) ~ harga,
is.na(harga) & !is.na(mean_pk) ~ mean_pk,
is.na(harga) & is.na(mean_pk) & !is.na(mean_pp) ~ mean_pp,
is.na(harga) & is.na(mean_pk) & is.na(mean_pp) & !is.na(mean_p) ~ mean_p,
TRUE ~ global_mean
))
sum(is.na(df_imp$harga_imputed)) # harus 0## [1] 0Agregasi per Provinsi
prov_prod <- df_imp %>%
group_by(nama_prov, produk) %>%
summarise(mean_harga = mean(harga_imputed, na.rm = TRUE)) %>%
ungroup() %>%
spread(key = produk, value = mean_harga)## `summarise()` has grouped output by 'nama_prov'. You can override using the
## `.groups` argument.# Jika beberapa provinsi masih NA, isi dengan rata-rata baris/kolom sebagai fallback
prov_prod_mat <- prov_prod %>%
column_to_rownames(var = "nama_prov") %>%
as.matrix()
# fallback: isi dengan rata-rata provinsi
row_means <- rowMeans(prov_prod_mat, na.rm = TRUE)
col_means <- colMeans(prov_prod_mat, na.rm = TRUE)
for(i in seq_len(nrow(prov_prod_mat))){
for(j in seq_len(ncol(prov_prod_mat))){
if(is.na(prov_prod_mat[i,j])){
prov_prod_mat[i,j] <- if(!is.na(row_means[i])) row_means[i] else col_means[j]
}
}
}df_agregasi <- df %>%
group_by(nama_prov, produk) %>%
summarise(rata_harga = mean(harga, na.rm = TRUE), .groups = "drop")## Warning: There were 510 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `rata_harga = mean(harga, na.rm = TRUE)`.
## ℹ In group 1: `nama_prov = "Aceh"` and `produk = "Bawang Merah"`.
## Caused by warning in `mean.default()`:
## ! argument is not numeric or logical: returning NA
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 509 remaining warnings.head(df_agregasi, 10)## # A tibble: 10 × 3
## nama_prov produk rata_harga
## <chr> <chr> <dbl>
## 1 Aceh Bawang Merah NA
## 2 Aceh Bawang Putih (Bonggol) NA
## 3 Aceh Beras Medium NA
## 4 Aceh Beras Premium NA
## 5 Aceh Cabai Merah Besar NA
## 6 Aceh Cabai Merah Keriting NA
## 7 Aceh Cabai Rawit Merah NA
## 8 Aceh Daging Ayam Ras NA
## 9 Aceh Daging Sapi Murni NA
## 10 Aceh Gula Konsumsi NAStandarisasi dan PCA
library(stats)
X <- scale(prov_prod_mat) # standar
if(ncol(X) > 1){
pca_res <- prcomp(X, center = FALSE, scale. = FALSE)
pca_scores <- pca_res$x[,1:2]
plot(pca_scores, main="PCA provinsi (PC1 vs PC2)")
} else {
# buat 2D pseudo untuk plotting (1 feature => use value and 0 for second)
pca_scores <- cbind(X[,1], rep(0, nrow(X)))
}
rownames(pca_scores) <- rownames(prov_prod_mat)Validasi Silhouette
sil_scores <- map_dbl(2:6, function(k){
km <- kmeans(X, centers = k, nstart = 25)
ss <- silhouette(km$cluster, dist(X))
mean(ss[,3], na.rm = TRUE)
})
sil_df <- data.frame(k = 2:6, silhouette = sil_scores)
print(sil_df)## k silhouette
## 1 2 0.6486772
## 2 3 0.7904051
## 3 4 0.6439106
## 4 5 0.5137600
## 5 6 0.4399051best_k <- sil_df$k[which.max(sil_df$silhouette)]
cat("Best k:", best_k, "with silhouette", max(sil_df$silhouette), "\n")## Best k: 3 with silhouette 0.7904051Silhouette score tertinggi diperoleh pada k = 3 dengan nilai 0.7904. Nilai ini menunjukkan struktur klaster yang sangat kuat sehingga penggunaan k = 3 dipertahankan untuk analisis selanjutnya.
Klasterisasi: K-Means & K-Medoids (PAM)
# KMeans
set.seed(123)
km_final <- kmeans(X, centers = best_k, nstart = 50)
km_labels <- km_final$cluster
# KMedoids (PAM)
pam_res <- pam(X, k = best_k, diss = FALSE)
pam_labels <- pam_res$clustering
# Tambahkan hasil ke data frame
res_df <- data.frame(
prov = rownames(prov_prod_mat),
km_cluster = km_labels,
pam_cluster = pam_labels,
pc1 = pca_scores[,1],
pc2 = pca_scores[,2],
stringsAsFactors = FALSE
)prov_names <- rownames(prov_prod_mat)
summary_df <- data.frame(prov = prov_names,
km_cluster = km_labels,
harga_daging = prov_prod_mat[, "Daging Sapi Murni"])
cluster_summary <- summary_df %>%
group_by(km_cluster) %>%
summarise(
n_provinsi = n(),
median_harga = median(harga_daging, na.rm = TRUE),
mean_harga = mean(harga_daging, na.rm = TRUE),
sd_harga = sd(harga_daging, na.rm = TRUE),
provinsi_list = paste(prov, collapse = ", "),
.groups = "drop"
) %>% arrange(km_cluster)
print(cluster_summary)## # A tibble: 3 × 6
## km_cluster n_provinsi median_harga mean_harga sd_harga provinsi_list
## <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 1 2 153936. 153936. 955. Aceh, Kepulauan Riau
## 2 2 5 134816. 135342. 1212. Bengkulu, Jambi, Lampu…
## 3 3 3 144190. 144320. 1900. Kepulauan Bangka Belit…# Heatmap (provinsi x fitur)
library(reshape2)
prov_df_long <- melt(prov_prod_mat)
colnames(prov_df_long) <- c("nama_prov", "produk", "mean_harga")
ggplot(prov_df_long, aes(x = produk, y = nama_prov, fill = mean_harga)) +
geom_tile() + scale_fill_viridis_c() +
labs(title = "Heatmap Harga (Provinsi x Produk)") +
theme_minimal() + theme(axis.text.x = element_text(angle=45,hjust=1))
# PCA scatter colored by KMeans
ggplot(res_df, aes(x = pc1, y = pc2, color = factor(km_cluster), label = prov)) +
geom_point(size=3) + geom_text(aes(label=prov), vjust=-1, size=3) +
labs(title = paste("PCA - KMeans (k=",best_k,")"), color="Cluster") + theme_minimal()
# PCA scatter colored by PAM
ggplot(res_df, aes(x = pc1, y = pc2, color = factor(pam_cluster), label = prov)) +
geom_point(size=3) + geom_text(aes(label=prov), vjust=-1, size=3) +
labs(title = paste("PCA - PAM (k=",best_k,")"), color="Cluster") + theme_minimal()
# Boxplot untuk Daging Sapi Murni per KMeans cluster
prod_focus <- "Daging Sapi Murni"
if(prod_focus %in% colnames(prov_prod_mat)){
box_df <- data.frame(prov = rownames(prov_prod_mat),
harga = prov_prod_mat[,prod_focus],
cluster = factor(km_labels))
ggplot(box_df, aes(x=cluster, y=harga)) + geom_boxplot() +
labs(title = paste("Distribusi harga",prod_focus,"per KMeans cluster"), x="Cluster", y="Harga")
}
- Klaster 1 (n = 2 provinsi): menunjukkan median
harga
Rp153.936/Kg, meliputi Provinsi Aceh dan Kepulauan Riau. Kedua provinsi ini cenderung merupakan provinsi pesisir sehingga harga relatif lebih tinggi kemungkinan karna biaya transportasi. - Klaster 2 (n = 5 provinsi): menunjukkan median
Rp134.816/Kgyang meliputi Provinsi Bengkulu, Jambi, Lampung, Sumatera Selatan, dan Sumatera Utara. Mayoritas provinsi dekat pasar besar/ibukota provinsi sehingga persaingan pasokan dan harga lebih rendah. - Klaster 3 (n = 3 provinsi): pada klaster ini median
harganya adalah
Rp144.190, meliputi Provinsi Kepulauan Bangka Belitung, Riau, dan Sumatera Barat. Ketiga provinsi ini harga daging sapi menengah karena permintaan yang stabil sehingga variasi harga antardaerah cenderung menengah.
Keterbatasan
Analisis ini memiliki beberapa keterbatasan. Pertama, imputasi mengisi nilai yang hilang dengan rata-rata wilayah sehingga dapat menurunkan varians asli dan menyamakan beberapa observasi yang seharusnya berbeda; hal ini dapat melemahkan kemampuan deteksi outlier. Kedua, model klaster tidak mempertimbangkan faktor penentu harga lain seperti biaya logistik (jarak ke pelabuhan/pasar), produksi lokal, tingkat konsumsi, atau kebijakan daerah — semua faktor yang berpengaruh pada harga daging.
Saran
- Melakukan clustering time series dengan metode seperti K-Medoids Fuzzy.
- Menambah variabel pendukung seperti jarak ke pelabuhan/pasar, produksi lokal, tingkat urbanisasi, dsb untuk menjelaskan perbedaan harga.