Analisis Klastering Dinamis Harga Daging Ayam Ras Menggunakan Dynamic Time Warping (DTW) di Provinsi Jawa Barat Tahun 2022–2025

Library

library(readxl)
library(tidyverse)
library(lubridate)
library(dplyr)
library(tidyr)
library(dtwclust)
library(ggplot2)
library(viridis)
library(sf)
library(plotly)

Import Data

path <- "D:/Github/STA1342-TPG/data_project_tpg.xlsx"
raw <- read_excel(path)
raw

## # A tibble: 27 × 47
##    region        `2022-1` `2022-2` `2022-3` `2022-4` `2022-5` `2022-6` `2022-7`
##    <chr>            <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
##  1 Bandung          35100    35111    32517    38300    38733    38577    36966
##  2 Bandung Barat    36484    34583    34409    36633    37422    37684    36796
##  3 Bekasi           39200    40000    40000    43143    47708    43042    43300
##  4 Bogor            37897    36643    35483    39333    39250    38704    37846
##  5 Ciamis           39800    37000    37933    43367    42500    37593    39774
##  6 Cianjur          36602    32958    34365    38589    36957    37448    36936
##  7 Cirebon          33871    30143    31000    34933    35400    34000    33433
##  8 Garut            35742    31762    32933    36044    34642    35893    34655
##  9 Indramayu        36133    31893    33645    36833    35226    38786    36742
## 10 Karawang         38536    36571    33567    31167    34714    35069    35633
## # ℹ 17 more rows
## # ℹ 39 more variables: `2022-8` <dbl>, `2022-9` <dbl>, `2022-10` <dbl>,
## #   `2022-11` <dbl>, `2022-12` <dbl>, `2023-1` <dbl>, `2023-2` <dbl>,
## #   `2023-3` <dbl>, `2023-4` <dbl>, `2023-5` <dbl>, `2023-6` <dbl>,
## #   `2023-7` <dbl>, `2023-8` <dbl>, `2023-9` <dbl>, `2023-10` <dbl>,
## #   `2023-11` <dbl>, `2023-12` <dbl>, `2024-1` <dbl>, `2024-2` <dbl>,
## #   `2024-3` <dbl>, `2024-4` <dbl>, `2024-5` <dbl>, `2024-6` <dbl>, …

peta <- st_read("D:\\Github\\idn_adm_bps_20200401_shp\\idn_admbnda_adm2_bps_20200401.shp")

## Reading layer `idn_admbnda_adm2_bps_20200401' from data source 
##   `D:\Github\idn_adm_bps_20200401_shp\idn_admbnda_adm2_bps_20200401.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 522 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84

peta <- peta[peta$ADM1_EN %in% c(
   "Jawa Barat"
), ]
buang <- c("Waduk Cirata")
peta <- subset(peta, !(ADM2_EN %in% buang))

Clustering

df <- raw %>%
  pivot_longer(
    cols = -region,
    names_to = "date",
    values_to = "price"
  )

df <- df %>%
  mutate(
    date = paste0(date, "-01"),
    date = as.Date(date),
    price = as.numeric(price)
  ) %>%
  arrange(region, date)

tsmat <- df %>%
  select(region, date, price) %>%
  pivot_wider(
    names_from = date,
    values_from = price
  ) %>%
  select(-region) %>%
  as.matrix()

library(dtwclust)

ks <- 2:8     # pilihan jumlah klaster
models <- list()

for (k in ks) {
  models[[as.character(k)]] <- tsclust(
    tsmat,
    type = "partitional",
    k = k,
    distance = "dtw_basic",
    centroid = "dba",
    seed = 123,
    trace = FALSE
  )
}

cvi_scores <- sapply(models, cvi)
cvi_scores

##                 2          3          4         5          6          7
## Sil     0.4720381  0.3281664  0.2390900 0.1188992 0.06699933 0.07613047
## SF      0.0000000  0.0000000  0.0000000 0.0000000 0.00000000 0.00000000
## CH     25.3524711 14.0923058 10.0153633 7.6363084 6.33400111 5.48192063
## DB      0.8639268  0.9122013  1.3513484 1.4384670 1.46606334 1.21242458
## DBstar  0.8639268  0.9322269  1.4936890 1.6633236 1.75734384 1.38399152
## D       0.3059116  0.4190948  0.3933283 0.3433045 0.34990836 0.39879000
## COP     0.3532514  0.2367529  0.2113031 0.2024202 0.19927646 0.17839297
##                 8
## Sil    0.03455885
## SF     0.00000000
## CH     4.74026976
## DB     1.10933986
## DBstar 1.43750279
## D      0.34111502
## COP    0.17632250

cluster_result <- tsclust(
  tsmat,
  type = "partitional",
  k = 2,
  distance = "dtw_basic",
  centroid = "dba",
  seed = 123,
  trace = TRUE
)

## Iteration 1: Changes / Distsum = 27 / 2103169
## Iteration 2: Changes / Distsum =  2 / 1594367
## Iteration 3: Changes / Distsum =  0 / 1549371
## 
##  Elapsed time is 0.06 seconds.

hasil <- data.frame(
  region = unique(df$region),
  cluster = cluster_result@cluster
)

hasil_wrap <- hasil %>%
  group_by(cluster) %>%
  summarise(region = list(region)) %>%
  unnest_longer(region)

hasil_wrap

## # A tibble: 27 × 2
##    cluster region       
##      <int> <chr>        
##  1       1 Bandung      
##  2       1 Bandung Barat
##  3       1 Bogor        
##  4       1 Cianjur      
##  5       1 Cirebon      
##  6       1 Garut        
##  7       1 Indramayu    
##  8       1 Karawang     
##  9       1 Kota Bandung 
## 10       1 Kota Banjar  
## # ℹ 17 more rows

Visualization

hc <- tsclust(
  tsmat,
  type = "hierarchical",
  distance = "dtw_basic",
  control = hierarchical_control(method = "average")
)

# Plot dendrogram
plot(hc)

peta_cluster <- peta %>%
  left_join(hasil, by = c("ADM2_EN" = "region"))

ggplot(peta_cluster) +
  geom_sf(aes(fill = factor(cluster)), color = "white", size = 0.2) +
  scale_fill_viridis_d(option = "plasma", begin = 0.2, end = 0.95,
                       name = "Cluster") +
  theme_minimal() +
  labs(title = "Choropleth Map Cluster")

df_plot <- df %>%
  left_join(hasil, by = "region")

ggplot(df_plot, aes(x = factor(cluster), y = price, fill = factor(cluster))) +
  geom_boxplot(alpha = 0.8) +
  scale_fill_viridis_d(
    option = "plasma",
    begin = 0.2,
    end = 0.95,
    name = "Cluster"
  ) +
  labs(
    title = "Distribusi Harga per Cluster",
    x = "Cluster",
    y = "Harga"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold")
  )

ggplot(df_plot, aes(x = date, y = price, group = region, color = factor(cluster))) +
  geom_line(alpha = 0.5, linewidth = 0.6) +
  scale_color_viridis_d(
    option = "plasma",   # palette viridis paling cerah
    begin = 0.2,         # mulai dari warna terang
    end = 0.95,          # hindari warna terlalu gelap
    name = "Cluster"
  ) +
  labs(
    title = "Pola Harga Daging Ayam per Wilayah",
    x = "Waktu",
    y = "Harga"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "right",
    legend.key.size = unit(0.7, "cm")
  )

cluster_mean <- df_plot %>%
  group_by(cluster, date) %>%
  summarise(mean_price = mean(price), .groups = "drop")

p <- ggplot(cluster_mean, aes(x = date, y = mean_price, color = factor(cluster))) +
  geom_line(linewidth = 1.1) +
  scale_color_viridis_d(option = "plasma", begin = 0.2, end = 0.95) +
  theme_minimal(base_size = 15)

ggplotly(p)

Summary

Berdasarkan hasil analis dan plot visualisasi, harga daging ayam ras di provinsi jawa barat pada tahun 2022-2025 beberapa dipengaruhi faktor geografis, hari perayaan, faktor eksternal (kebijakan pemerintah, musibah/krisis produksi, kenaikan harga pakan ayam, dll)