setwd("C:/Users/sambe/Downloads/TUGAS 1 STAT LINGKUNGAN")
getwd()

file_data <- "Rata-Rata_Volume Sampah Sungai di Provinsi DKI Jakarta Menurut Kota Administrasi.xls"

tabel_data <- read_html(file_data) |>
  html_table(fill = TRUE)

data_raw <- tabel_data[[1]]

data <- data_raw |>
  clean_names() |>
  mutate(
    periode_data = as.character(periode_data),
    tanggal = ymd(paste0(periode_data, "01")),
    wilayah = str_squish(wilayah),
    kecamatan = str_squish(kecamatan),
    volume_sampah = as.numeric(volume_sampah),
    total_sampah_residu = as.numeric(total_sampah_residu),
    rata_rata = as.numeric(rata_rata)
  ) |>
  filter(!is.na(rata_rata))

glimpse(data)

data |>
  summarise(
    jumlah_observasi = n(),
    jumlah_kecamatan = n_distinct(kecamatan),
    jumlah_kota = n_distinct(wilayah),
    periode_awal = min(tanggal, na.rm = TRUE),
    periode_akhir = max(tanggal, na.rm = TRUE)
  ) |>
  kable(caption = "Ringkasan Struktur Data") |>
  kable_styling(full_width = FALSE)

data |>
  count(wilayah, sort = TRUE) |>
  kable(
    col.names = c("Kota Administrasi", "Jumlah Observasi"),
    caption = "Jumlah Observasi Berdasarkan Kota Administrasi"
  ) |>
  kable_styling(full_width = FALSE)

statistik_deskriptif <- data |>
  summarise(
    N = n(),
    Mean = mean(rata_rata, na.rm = TRUE),
    Median = median(rata_rata, na.rm = TRUE),
    SD = sd(rata_rata, na.rm = TRUE),
    Varians = var(rata_rata, na.rm = TRUE),
    Minimum = min(rata_rata, na.rm = TRUE),
    Q1 = quantile(rata_rata, 0.25, na.rm = TRUE),
    Q3 = quantile(rata_rata, 0.75, na.rm = TRUE),
    IQR = IQR(rata_rata, na.rm = TRUE),
    Maksimum = max(rata_rata, na.rm = TRUE)
  )

statistik_deskriptif |>
  mutate(across(where(is.numeric), ~ round(.x, 2))) |>
  kable(caption = "Statistik Deskriptif Volume Sampah Sungai") |>
  kable_styling(full_width = FALSE)

statistik_kota <- data |>
  group_by(wilayah) |>
  summarise(
    N = n(),
    Mean = mean(rata_rata, na.rm = TRUE),
    Median = median(rata_rata, na.rm = TRUE),
    SD = sd(rata_rata, na.rm = TRUE),
    Minimum = min(rata_rata, na.rm = TRUE),
    Q1 = quantile(rata_rata, 0.25, na.rm = TRUE),
    Q3 = quantile(rata_rata, 0.75, na.rm = TRUE),
    Maksimum = max(rata_rata, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(across(where(is.numeric), ~ round(.x, 2)))

statistik_kota |>
  kable(
    col.names = c("Kota Administrasi", "N", "Mean", "Median", "SD",
                  "Minimum", "Q1", "Q3", "Maksimum"),
    caption = "Statistik Deskriptif Berdasarkan Kota Administrasi"
  ) |>
  kable_styling(full_width = FALSE)

ggplot(data, aes(x = rata_rata)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_vline(
    aes(xintercept = mean(rata_rata, na.rm = TRUE)),
    linetype = "dashed", linewidth = 1
  ) +
  labs(
    title = "Distribusi Rata-Rata Volume Sampah Sungai",
    subtitle = "DKI Jakarta, 2024–2026",
    x = "Rata-rata Volume Sampah",
    y = "Frekuensi"
  ) +
  theme_minimal()

ggplot(data, aes(y = rata_rata)) +
  geom_boxplot(fill = "lightblue", outlier.alpha = 0.5) +
  labs(
    title = "Boxplot Rata-Rata Volume Sampah Sungai",
    x = NULL,
    y = "Rata-rata Volume Sampah"
  ) +
  theme_minimal()

ggplot(
  data,
  aes(
    x = reorder(wilayah, rata_rata, FUN = median),
    y = rata_rata
  )
) +
  geom_boxplot(aes(group = wilayah), fill = "lightblue", outlier.alpha = 0.4) +
  labs(
    title = "Distribusi Volume Sampah Berdasarkan Kota Administrasi",
    x = "Kota Administrasi",
    y = "Rata-rata Volume Sampah"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 35, hjust = 1))

data_bulanan <- data |>
  group_by(tanggal) |>
  summarise(
    rata_rata_volume = mean(rata_rata, na.rm = TRUE),
    median_volume = median(rata_rata, na.rm = TRUE),
    minimum = min(rata_rata, na.rm = TRUE),
    maksimum = max(rata_rata, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(data_bulanan, aes(x = tanggal, y = rata_rata_volume)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.8) +
  labs(
    title = "Perubahan Rata-Rata Volume Sampah Sungai",
    x = "Periode",
    y = "Rata-rata Volume Sampah"
  ) +
  scale_x_date(date_breaks = "3 months", date_labels = "%b\n%Y") +
  theme_minimal()

data_kota_bulanan <- data |>
  group_by(wilayah, tanggal) |>
  summarise(
    rata_rata_volume = mean(rata_rata, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(
  data_kota_bulanan,
  aes(x = tanggal, y = rata_rata_volume, group = wilayah)
) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 1.3) +
  facet_wrap(~ wilayah, scales = "free_y") +
  labs(
    title = "Pola Temporal Volume Sampah Berdasarkan Kota Administrasi",
    x = "Periode",
    y = "Rata-rata Volume Sampah"
  ) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b\n%Y") +
  theme_minimal()

Q1 <- quantile(data$rata_rata, 0.25, na.rm = TRUE)
Q3 <- quantile(data$rata_rata, 0.75, na.rm = TRUE)
IQR_value <- IQR(data$rata_rata, na.rm = TRUE)

batas_bawah <- Q1 - 1.5 * IQR_value
batas_atas <- Q3 + 1.5 * IQR_value

tibble(
  Q1 = Q1,
  Q3 = Q3,
  IQR = IQR_value,
  Batas_Bawah = batas_bawah,
  Batas_Atas = batas_atas
) |>
  mutate(across(where(is.numeric), ~ round(.x, 2))) |>
  kable(caption = "Batas Outlier Berdasarkan Metode IQR") |>
  kable_styling(full_width = FALSE)

data_outlier <- data |>
  mutate(
    status_outlier = case_when(
      rata_rata < batas_bawah ~ "Outlier bawah",
      rata_rata > batas_atas ~ "Outlier atas",
      TRUE ~ "Normal"
    )
  )

data_outlier |>
  count(status_outlier) |>
  kable(
    col.names = c("Status", "Jumlah"),
    caption = "Jumlah Observasi Berdasarkan Status Outlier"
  ) |>
  kable_styling(full_width = FALSE)

data_outlier |>
  filter(status_outlier != "Normal") |>
  select(periode_data, tanggal, wilayah, kecamatan, rata_rata, status_outlier) |>
  arrange(desc(rata_rata)) |>
  head(20) |>
  mutate(rata_rata = round(rata_rata, 2)) |>
  kable(
    col.names = c("Periode", "Tanggal", "Kota Administrasi",
                  "Kecamatan", "Rata-rata", "Status"),
    caption = "20 Observasi Outlier dengan Nilai Tertinggi"
  ) |>
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover")
  )

ggplot(data_outlier, aes(x = tanggal, y = rata_rata, color = status_outlier)) +
  geom_point(alpha = 0.65) +
  geom_hline(yintercept = batas_atas, linetype = "dashed", linewidth = 0.8) +
  labs(
    title = "Identifikasi Outlier Volume Sampah Sungai",
    x = "Periode",
    y = "Rata-rata Volume Sampah",
    color = "Status"
  ) +
  theme_minimal()

acf(
  data_bulanan$rata_rata_volume,
  main = "Autocorrelation Function Rata-Rata Volume Sampah"
)

acf_result <- acf(
  data_bulanan$rata_rata_volume,
  plot = FALSE
)

tibble(
  lag = as.numeric(acf_result$lag),
  autocorrelation = as.numeric(acf_result$acf)
) |>
  filter(lag > 0) |>
  head(12) |>
  mutate(
    lag = round(lag, 0),
    autocorrelation = round(autocorrelation, 3)
  ) |>
  kable(caption = "Nilai Autokorelasi Temporal pada Beberapa Lag") |>
  kable_styling(full_width = FALSE)

folder_kota <- c(
  "JAKARTA PUSAT",
  "JAKARTA SELATAN",
  "JAKARTA BARAT",
  "JAKARTA TIMUR",
  "JAKARTA UTARA"
)

file_shp <- file.path(
  folder_kota,
  "ADMINISTRASIKECAMATAN_AR_25K.shp"
)

tibble(
  kota = folder_kota,
  file_shp = file_shp,
  tersedia = file.exists(file_shp)
) |>
  kable(caption = "Daftar Shapefile Batas Kecamatan") |>
  kable_styling(full_width = FALSE)

shp_list <- map2(
  file_shp,
  folder_kota,
  ~ st_read(.x, quiet = TRUE) |>
    st_make_valid() |>
    st_transform(4326) |>
    mutate(kota_shp = .y)
)

shp_jakarta <- bind_rows(shp_list)

tibble(
  nomor = seq_along(names(shp_jakarta)),
  field = names(shp_jakarta)
) |>
  kable(caption = "Field pada Shapefile Administrasi Kecamatan") |>
  kable_styling(full_width = FALSE)

kecamatan_data <- data |>
  mutate(kecamatan_clean = str_to_upper(str_squish(kecamatan))) |>
  distinct(kecamatan_clean) |>
  pull(kecamatan_clean)

candidate_fields <- names(shp_jakarta)[
  !names(shp_jakarta) %in% c("geometry", "kota_shp")
]

field_score <- map_dfr(
  candidate_fields,
  function(f) {
    nilai_shp <- shp_jakarta[[f]]

    if (!is.character(nilai_shp) && !is.factor(nilai_shp)) {
      return(tibble(field = f, matched = 0, total_unique = 0))
    }

    nilai_shp <- str_to_upper(str_squish(as.character(nilai_shp)))

    tibble(
      field = f,
      matched = sum(nilai_shp %in% kecamatan_data, na.rm = TRUE),
      total_unique = n_distinct(nilai_shp)
    )
  }
) |>
  arrange(desc(matched), desc(total_unique))

field_score |>
  head(15) |>
  kable(
    caption = "Skor Kecocokan Field Shapefile dengan Nama Kecamatan"
  ) |>
  kable_styling(full_width = FALSE)

field_kecamatan <- field_score |>
  filter(matched == max(matched, na.rm = TRUE)) |>
  slice(1) |>
  pull(field)

field_kecamatan

shp_jakarta <- shp_jakarta |>
  mutate(
    kecamatan = str_to_upper(
      str_squish(as.character(.data[[field_kecamatan]]))
    )
  )

data_spasial <- data |>
  mutate(
    kecamatan = str_to_upper(str_squish(kecamatan))
  )

shp_jakarta |>
  st_drop_geometry() |>
  distinct(kecamatan, kota_shp) |>
  mutate(cocok = kecamatan %in% data_spasial$kecamatan) |>
  count(cocok) |>
  kable(
    col.names = c("Nama Kecamatan Cocok", "Jumlah"),
    caption = "Pemeriksaan Kecocokan Nama Kecamatan"
  ) |>
  kable_styling(full_width = FALSE)

jakarta_map |>
  st_drop_geometry() |>
  summarise(
    total_polygon = n(),
    polygon_dengan_data = sum(!is.na(rata_rata)),
    polygon_tanpa_data = sum(is.na(rata_rata))
  ) |>
  kable(caption = "Pemeriksaan Penggabungan Data Statistik dan Spasial") |>
  kable_styling(full_width = FALSE)

ggplot(
  jakarta_map |> filter(!is.na(rata_rata)),
  aes(fill = rata_rata)
) +
  geom_sf(color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "C",
    name = "Rata-rata\nvolume"
  ) +
  facet_wrap(~ kota_shp) +
  labs(
    title = "Distribusi Spasial Rata-Rata Volume Sampah Sungai",
    subtitle = "Rata-rata seluruh periode pada tingkat kecamatan",
    x = NULL, y = NULL
  ) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    panel.grid = element_blank(),
    strip.text = element_text(face = "bold")
  )

data_spasial_kecamatan |>
  arrange(desc(rata_rata)) |>
  mutate(peringkat = row_number()) |>
  select(peringkat, wilayah, kecamatan, rata_rata, median, maksimum) |>
  head(15) |>
  mutate(across(where(is.numeric), ~ round(.x, 2))) |>
  kable(
    col.names = c("Peringkat", "Kota Administrasi", "Kecamatan",
                  "Mean", "Median", "Maksimum"),
    caption = "15 Kecamatan dengan Rata-Rata Volume Sampah Tertinggi"
  ) |>
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover")
  )

moran_map <- jakarta_map |>
  filter(!is.na(rata_rata), !is.na(kota_shp)) |>
  arrange(kota_shp, kecamatan)

nb_jakarta <- poly2nb(
  moran_map,
  queen = TRUE
)

lw_jakarta <- nb2listw(
  nb_jakarta,
  style = "W",
  zero.policy = TRUE
)

tibble(
  kecamatan = moran_map$kecamatan,
  jumlah_tetangga = card(nb_jakarta)
) |>
  arrange(desc(jumlah_tetangga)) |>
  kable(caption = "Jumlah Tetangga Setiap Kecamatan") |>
  kable_styling(full_width = FALSE)

moran_result <- moran.test(
  moran_map$rata_rata,
  lw_jakarta,
  zero.policy = TRUE
)

moran_result

tibble(
  Moran_I = as.numeric(moran_result$estimate[1]),
  Expected_I = as.numeric(moran_result$estimate[2]),
  Variance = as.numeric(moran_result$estimate[3]),
  P_value = moran_result$p.value
) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "Ringkasan Global Moran's I") |>
  kable_styling(full_width = FALSE)

nilai_rata2 <- mean(moran_map$rata_rata, na.rm = TRUE)

moran_map <- moran_map |>
  mutate(
    lag_rata_rata = lag.listw(
      lw_jakarta,
      rata_rata,
      zero.policy = TRUE
    ),
    cluster = case_when(
      rata_rata >= nilai_rata2 &
        lag_rata_rata >= nilai_rata2 &
        local_p < 0.05 ~ "High-High",

      rata_rata < nilai_rata2 &
        lag_rata_rata < nilai_rata2 &
        local_p < 0.05 ~ "Low-Low",

      rata_rata >= nilai_rata2 &
        lag_rata_rata < nilai_rata2 &
        local_p < 0.05 ~ "High-Low",

      rata_rata < nilai_rata2 &
        lag_rata_rata >= nilai_rata2 &
        local_p < 0.05 ~ "Low-High",

      TRUE ~ "Tidak signifikan"
    )
  )

moran_map |>
  st_drop_geometry() |>
  count(cluster) |>
  arrange(desc(n)) |>
  kable(
    col.names = c("Kategori Cluster", "Jumlah Kecamatan"),
    caption = "Klasifikasi Local Spatial Autocorrelation"
  ) |>
  kable_styling(full_width = FALSE)