1 Pendahuluan

1.1 Tentang Studi Kasus

Dokumen ini menerapkan seluruh teknik eksplorasi data lingkungan yang dibahas pada Pertemuan 2 — statistik deskriptif, visualisasi, deteksi outlier, dan autokorelasi (temporal & spasial) — pada studi kasus kualitas udara di Provinsi Banten.

Dataset yang digunakan adalah data simulasi konsentrasi PM2.5 (µg/m³) rata-rata bulanan dan curah hujan (mm) dari 6 stasiun pemantauan di Provinsi Banten (Cilegon, Kota Serang, Kota Tangerang, Tangerang Selatan, Pandeglang, dan Lebak) sepanjang Januari–Desember 2025 (72 baris = 6 stasiun × 12 bulan), sebagaimana dijelaskan pada slide materi Pertemuan 2. Download data disini (link data)

1.2 Rumusan Pertanyaan Analisis

  1. Bagaimana karakteristik statistik deskriptif konsentrasi PM2.5 secara keseluruhan dan per stasiun?
  2. Bagaimana pola distribusi, tren bulanan, dan sebaran spasial PM2.5 dapat digambarkan melalui visualisasi?
  3. Apakah terdapat episode konsentrasi PM2.5 yang tergolong outlier? Di stasiun dan bulan apa?
  4. Apakah terdapat indikasi autokorelasi temporal (antarbulan) dan autokorelasi spasial (antarstasiun berdekatan) pada data PM2.5?

Keempat pertanyaan ini dijawab bertahap melalui delapan langkah analisis di bawah ini.

2 Persiapan

2.1 Memuat Package

library(readxl)   # membaca file Excel
library(dplyr)    # manipulasi data
library(ggplot2)  # visualisasi
library(sf)       # objek spasial
library(tmap)     # peta tematik
library(spdep)    # autokorelasi spasial (Moran's I)

2.2 Langkah 1 — Import dan Struktur Data

# Membaca data kualitas udara Provinsi Banten
# (sesuaikan path bila file tidak berada di working directory yang sama)
setwd("C:/Users/mifta/OneDrive - untirta.ac.id/UNTIRTA/STATISTIKA/Materi MK Semester Ganjil/Statistika Lingkungan")
data_pm25 <- read_excel("2-Data Simulasi - Kualitas Udara Banten 2025.xlsx")

# Melihat struktur dan ukuran data
str(data_pm25)
## tibble [72 × 10] (S3: tbl_df/tbl/data.frame)
##  $ stasiun         : chr [1:72] "Cilegon" "Cilegon" "Cilegon" "Cilegon" ...
##  $ kabupaten_kota  : chr [1:72] "Kota Cilegon" "Kota Cilegon" "Kota Cilegon" "Kota Cilegon" ...
##  $ kategori_kawasan: chr [1:72] "Kawasan Industri" "Kawasan Industri" "Kawasan Industri" "Kawasan Industri" ...
##  $ bulan           : num [1:72] 1 2 3 4 5 6 7 8 9 10 ...
##  $ nama_bulan      : chr [1:72] "Januari" "Februari" "Maret" "April" ...
##  $ tahun           : num [1:72] 2025 2025 2025 2025 2025 ...
##  $ pm25            : num [1:72] 19 20.6 18.6 22.9 28.3 29.3 33.3 33.5 96.4 28.5 ...
##  $ curah_hujan     : num [1:72] 288 256 243 165 112 ...
##  $ lat             : num [1:72] -5.99 -5.99 -5.99 -5.99 -5.99 ...
##  $ long            : num [1:72] 106 106 106 106 106 ...
head(data_pm25, 6)
dim(data_pm25)
## [1] 72 10
# Memeriksa data hilang (missing value)
sum(is.na(data_pm25))
## [1] 0

Interpretasi

Dataset berisi 72 baris dan 10 kolom, sesuai rancangan 6 stasiun × 12 bulan. Kolom pm25 dan curah_hujan sudah bertipe numerik, stasiun bertipe karakter, dan tidak ditemukan nilai hilang (sum(is.na(...)) = 0) — artinya data siap dianalisis tanpa perlu tahap imputasi. Kolom lat/long akan dipakai pada Langkah 5 dan 8, dan kolom bulan (numerik 1–12) akan dipakai pada Langkah 4 dan 7.

3 Statistik Deskriptif

3.1 Langkah 2 — Statistik Deskriptif

# Statistik deskriptif keseluruhan
summary(data_pm25$pm25)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.60   18.52   24.95   28.14   33.80   96.40
sd(data_pm25$pm25)
## [1] 14.45932
IQR(data_pm25$pm25)
## [1] 15.275
# Statistik deskriptif per stasiun
rekap_stasiun <- data_pm25 %>%
  group_by(stasiun) %>%
  summarise(
    n      = n(),
    rerata = mean(pm25),
    median = median(pm25),
    sd     = sd(pm25),
    min    = min(pm25),
    maks   = max(pm25)
  ) %>%
  arrange(desc(rerata))

rekap_stasiun

Interpretasi

Secara keseluruhan, rata-rata PM2.5 adalah 28,1 µg/m³ dengan median 25,0 µg/m³ dan simpangan baku 14,5 µg/m³; rentang antarkuartil (IQR) sebesar 15,3. Selisih antara mean dan median ini adalah petunjuk awal bahwa distribusi condong ke kanan (right-skewed) — beberapa nilai yang sangat tinggi menarik rata-rata ke atas melebihi median.

Rekap per stasiun memperlihatkan gradien yang jelas:

  • Kota Tangerang (rerata 42,7) dan Tangerang Selatan (37,2) — dua stasiun perkotaan padat — memiliki rerata PM2.5 tertinggi.
  • Cilegon (kawasan industri) berada di tengah (31,1), namun simpangan bakunya jauh lebih besar (sd = 21,2, dibanding stasiun lain yang berkisar 3–13) — ini mengindikasikan adanya nilai ekstrem di stasiun tersebut (akan dikonfirmasi di Langkah 6).
  • Kota Serang (24,3) berada di tingkat menengah, sedangkan Pandeglang (18,3) dan Lebak (15,2) — dua stasiun pedesaan — konsisten menjadi yang terendah.

Pola ini sejalan dengan intuisi: wilayah perkotaan padat dan kawasan industri cenderung memiliki beban pencemaran udara yang lebih tinggi dibanding wilayah pedesaan.

4 Visualisasi Data

4.1 Langkah 3 — Histogram dan Boxplot

# Histogram distribusi PM2.5 keseluruhan
ggplot(data_pm25, aes(x = pm25)) +
  geom_histogram(binwidth = 5, fill = "#147D64", color = "white") +
  labs(title = "Distribusi Konsentrasi PM2.5 - Provinsi Banten 2025",
       x = "PM2.5 (µg/m³)", y = "Frekuensi") +
  theme_minimal()

# Boxplot PM2.5 per stasiun
ggplot(data_pm25, aes(x = reorder(stasiun, pm25, median), y = pm25)) +
  geom_boxplot(fill = "#FFC200") +
  labs(title = "Sebaran PM2.5 per Stasiun",
       x = "Stasiun", y = "PM2.5 (µg/m³)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 25, hjust = 1))

Interpretasi

Histogram memperlihatkan sebagian besar pengamatan terkonsentrasi pada rentang 15–40 µg/m³, dengan ekor memanjang ke kanan menuju nilai 80–100 µg/m³ — mengonfirmasi kecondongan (skewness) positif yang sudah terlihat dari selisih mean-median pada Langkah 2.

Boxplot per stasiun (diurutkan dari median terendah ke tertinggi) menunjukkan:

  • Median Kota Tangerang dan Tangerang Selatan berada jauh di atas stasiun lain, sejalan dengan hasil statistik deskriptif.
  • Kotak (IQR) Cilegon dan Kota Tangerang tampak lebih pendek dibanding rentang whisker-nya, dengan titik terpisah di atas whisker atas — inilah representasi visual dari kedua nilai ekstrem yang akan diidentifikasi secara formal pada Langkah 6.
  • Pandeglang dan Lebak memiliki kotak yang relatif rendah dan sempit, menandakan konsentrasi PM2.5 yang stabil dan rendah sepanjang tahun.

4.2 Langkah 4 — Diagram Deret Waktu

ggplot(data_pm25, aes(x = bulan, y = pm25, color = stasiun)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  scale_x_continuous(breaks = 1:12) +
  labs(title = "Tren Bulanan PM2.5 Tahun 2025",
       x = "Bulan", y = "PM2.5 (µg/m³)", color = "Stasiun") +
  theme_minimal()

# Pembanding: rata-rata PM2.5 vs curah hujan per bulan (semua stasiun)
rerata_bulanan <- data_pm25 %>%
  group_by(bulan) %>%
  summarise(rerata_pm25 = mean(pm25), rerata_hujan = mean(curah_hujan))

# Membuat scatter plot dengan garis tren
ggplot(rerata_bulanan, aes(x = rerata_pm25, y = rerata_hujan)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "lm", color = "red", se = TRUE) + # se=TRUE menampilkan area confidence interval
  labs(title = "Scatter Plot Korelasi rata-rata PM2.5 vs curah hujan per bulan",
       x = "rata-rata PM2.5 (µg/m³)",
       y = "curah hujan (mm)") +
  theme_minimal()

cor(rerata_bulanan$rerata_pm25, rerata_bulanan$rerata_hujan)
## [1] -0.8430532

Interpretasi

Grafik deret waktu memperlihatkan pola musiman yang jelas dan konsisten di seluruh stasiun: PM2.5 relatif rendah pada Januari–April (rerata ± 22–24 µg/m³), meningkat tajam memasuki musim kemarau (Mei–September), memuncak pada bulan September (rerata lintas-stasiun 43,4 µg/m³), lalu turun kembali begitu musim hujan tiba (Oktober–Desember).

Korelasi antara rata-rata bulanan PM2.5 dan curah hujan bernilai negatif kuat (sekitar -0,84) — mengonfirmasi secara kuantitatif bahwa bulan-bulan dengan curah hujan tinggi cenderung memiliki kualitas udara yang lebih baik (curah hujan “membersihkan” partikulat di udara), sementara musim kemarau yang kering berasosiasi dengan konsentrasi PM2.5 yang lebih tinggi — plot juga menunjukkan lonjakan tajam di Cilegon (September) dan Kota Tangerang (Agustus) yang menonjol dibanding pola musiman stasiun lain.

4.3 Langkah 5 — Peta Tematik

# Rata-rata PM2.5 per stasiun beserta koordinatnya
rerata_stasiun <- data_pm25 %>%
  group_by(kabupaten_kota, lat, long) %>%
  summarise(rerata_pm25 = mean(pm25), .groups = "drop")

rerata_stasiun
#Import SHP Banten
shp<-st_read("2-SHPBANTEN/Banten_ADMIN_BPS.shp")
## Reading layer `Banten_ADMIN_BPS' from data source 
##   `C:\Users\mifta\OneDrive - untirta.ac.id\UNTIRTA\STATISTIKA\Materi MK Semester Ganjil\Statistika Lingkungan\2-SHPBANTEN\Banten_ADMIN_BPS.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -7.016779 xmax: 106.7799 ymax: -5.807418
## Geodetic CRS:  WGS 84
#Gabungkan Data shapefile kabupaten/kota
titik_stasiun <- shp %>%
  left_join(rerata_stasiun, by = c("Kabupaten" = "kabupaten_kota"))

#Pemetaan Variabel rerata_pm25
library(mapview)
mapview(titik_stasiun, zcol = "rerata_pm25")
# Ubah menjadi objek spasial
titik_stasiun <- st_as_sf(rerata_stasiun, coords = c("long", "lat"), crs = 4326)

# Peta tematik: ukuran & warna titik merepresentasikan rerata PM2.5
tm_shape(titik_stasiun) +
  tm_symbols(size = "rerata_pm25", col = "rerata_pm25",
             palette = "YlOrRd", title.col = "Rerata PM2.5",
             title.size = "Rerata PM2.5") +
  tm_layout(title = "Sebaran Kualitas Udara di Provinsi Banten (Rerata 2025)")

Interpretasi

Peta tematik memperlihatkan pola spasial yang cukup rapi: kluster timur (Kota Tangerang dan Tangerang Selatan — wilayah urban padat yang berbatasan dengan DKI Jakarta) tampil dengan simbol besar dan warna gelap (PM2.5 tinggi), sementara kluster barat/selatan (Pandeglang dan Lebak — wilayah pedesaan) tampil dengan simbol kecil dan warna terang (PM2.5 rendah). Cilegon dan Kota Serang berada pada posisi menengah. Pola “yang berdekatan cenderung mirip” ini merupakan petunjuk visual awal adanya autokorelasi spasial positif, yang akan diuji secara formal pada Langkah 8.

5 Deteksi Outlier

5.1 Langkah 6 — Deteksi Outlier

# Metode IQR
Q1 <- quantile(data_pm25$pm25, 0.25)
Q3 <- quantile(data_pm25$pm25, 0.75)
IQR_val <- Q3 - Q1
batas_bawah <- Q1 - 1.5 * IQR_val
batas_atas  <- Q3 + 1.5 * IQR_val

cat("Batas bawah:", round(batas_bawah, 1), " | Batas atas:", round(batas_atas, 1), "\n")
## Batas bawah: -4.4  | Batas atas: 56.7
outlier_iqr <- data_pm25 %>% filter(pm25 < batas_bawah | pm25 > batas_atas)
outlier_iqr
# Metode Z-score
data_pm25$z_pm25 <- as.numeric(scale(data_pm25$pm25))
outlier_z <- data_pm25 %>% filter(abs(z_pm25) > 3)
outlier_z %>% select(stasiun, bulan, nama_bulan, pm25, z_pm25)

Interpretasi

Dengan aturan IQR, batas atas kewajaran adalah 56,7 µg/m³ (batas bawah bernilai negatif sehingga tidak relevan untuk data konsentrasi). Dua pengamatan melampaui batas ini dan kedua metode (IQR dan Z-score) sepakat menandai titik yang sama:

Stasiun Bulan PM2.5 Z-score
Cilegon September 96,4 +4,72
Kota Tangerang Agustus 79,8 +3,57

Kesepakatan dua metode ini memperkuat keyakinan bahwa kedua nilai memang menyimpang signifikan dari pola umum data, bukan sekadar batas ambang yang kebetulan terlampaui. Karena keduanya terjadi tepat pada puncak musim kemarau — periode yang secara konsisten sudah menunjukkan tren meningkat pada Langkah 4 — nilai ini lebih masuk akal diinterpretasikan sebagai episode polusi ekstrem yang nyata (misalnya kabut asap akibat kebakaran lahan/hutan, yang umum terjadi di Indonesia saat kemarau panjang) dibandingkan kesalahan pengukuran. Sesuai prinsip yang dibahas di kelas, kedua titik ini tidak boleh langsung dibuang — sebaliknya, keduanya adalah temuan penting yang layak ditelusuri lebih lanjut (mis. dikonfirmasi dengan data titik panas/hotspot kebakaran pada bulan yang sama).

6 Autokorelasi

6.1 Langkah 7 — Autokorelasi Temporal (ACF)

# Ekstrak deret waktu PM2.5 stasiun Cilegon
ts_cilegon <- data_pm25 %>% filter(stasiun == "Cilegon") %>% pull(pm25)
ts_cilegon <- ts(ts_cilegon, start = c(2025, 1), frequency = 12)

acf(ts_cilegon, main = "ACF PM2.5 - Stasiun Cilegon")

# Pembanding: stasiun tanpa nilai ekstrem (Pandeglang)
ts_pandeglang <- data_pm25 %>% filter(stasiun == "Pandeglang") %>% pull(pm25)
ts_pandeglang <- ts(ts_pandeglang, start = c(2025, 1), frequency = 12)

acf(ts_pandeglang, main = "ACF PM2.5 - Stasiun Pandeglang (pembanding)")

Interpretasi

ACF pada Cilegon menunjukkan autokorelasi lag-1 yang relatif lemah (± 0,10) dan tidak melewati garis signifikansi — sekilas seperti “tidak ada pola”. Namun ini perlu dibaca hati-hati: pada Langkah 6 kita sudah tahu Cilegon memiliki outlier ekstrem di bulan September. Satu nilai ekstrem pada deret waktu yang pendek (hanya 12 titik) dapat mendistorsi/menutupi pola autokorelasi musiman yang sebenarnya ada.

Sebagai pembanding, ACF Pandeglang — stasiun dengan pola musiman yang mulus tanpa nilai ekstrem — menunjukkan autokorelasi lag-1 yang kuat dan positif (± 0,73), jauh melewati garis signifikansi. Perbandingan ini adalah contoh nyata mengapa outlier perlu ditangani/diperhatikan sebelum menafsirkan autokorelasi: kesimpulan “tidak ada autokorelasi” pada Cilegon bisa jadi keliru apabila semata didasarkan pada ACF mentah tanpa mempertimbangkan pencilan yang sudah teridentifikasi sebelumnya.

(Sebagai latihan tambahan, coba jalankan acf() untuk keempat stasiun lainnya dan bandingkan nilai lag-1-nya.)

6.2 Langkah 8 — Autokorelasi Spasial (Moran’s I)

# Matriks pembobot spasial: k = 3 tetangga terdekat
koordinat <- cbind(rerata_stasiun$long, rerata_stasiun$lat)
knn3 <- knearneigh(koordinat, k = 3)
wij <- nb2listw(knn2nb(knn3), style = "W")

# Melihat daftar tetangga tiap stasiun
knn2nb(knn3)
## Neighbour list object:
## Number of regions: 6 
## Number of nonzero links: 18 
## Percentage nonzero weights: 50 
## Average number of links: 3 
## Non-symmetric neighbours list
# Uji Moran's I
moran.test(rerata_stasiun$rerata_pm25, listw = wij)
## 
##  Moran I test under randomisation
## 
## data:  rerata_stasiun$rerata_pm25  
## weights: wij    
## 
## Moran I statistic standard deviate = 2.8764, p-value = 0.002011
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.28125027       -0.20000000        0.02799318

Interpretasi

Dengan skema k = 3 tetangga terdekat, hasil simulasi ini menghasilkan Indeks Moran’s I sekitar +0,28 dengan nilai-p ≈ 0,002 (< 0,05) — menunjukkan autokorelasi spasial positif dan signifikan secara statistik: stasiun yang berdekatan secara geografis cenderung memiliki tingkat PM2.5 yang mirip, bukan acak.

Daftar tetangga (knn2nb) memperjelas struktur klusternya: Pandeglang dan Lebak konsisten saling bertetangga bersama Cilegon dan Kota Serang (kluster barat/selatan dengan PM2.5 rendah-menengah), sementara Kota Tangerang dan Tangerang Selatan saling bertetangga sebagai kluster timur ber-PM2.5 tinggi — namun keduanya juga menjadikan Cilegon dan Kota Serang sebagai tetangga terdekat ketiga. Cilegon dan Kota Serang inilah yang berperan sebagai penghubung antarkluster, yang sedikit melemahkan (namun tidak menghilangkan) kekuatan pengelompokan spasial secara keseluruhan. Pola ini konsisten dengan Hukum Pertama Geografi Tobler yang dibahas di kelas: kedekatan lokasi berasosiasi dengan kedekatan nilai.

7 Kesimpulan

Menjawab keempat pertanyaan analisis di awal:

  1. Statistik deskriptif menunjukkan gradien PM2.5 yang jelas: perkotaan padat (Tangerang, Tangerang Selatan) > kawasan industri (Cilegon) > perkotaan (Serang) > pedesaan (Pandeglang, Lebak); distribusi keseluruhan condong ke kanan.
  2. Visualisasi mengonfirmasi pola musiman (memuncak September, korelasi negatif dengan curah hujan) dan pola spasial (kluster timur tinggi, kluster barat/selatan rendah).
  3. Deteksi outlier menemukan dua episode ekstrem (Cilegon September; Kota Tangerang Agustus) yang konsisten dengan narasi kabut asap musim kemarau, dan perlu ditelaah lebih lanjut alih-alih dibuang begitu saja.
  4. Autokorelasi temporal terbukti ada namun bervariasi antarstasiun (dan dapat terselubung oleh outlier, seperti pada kasus Cilegon), sementara autokorelasi spasial terbukti positif dan signifikan (Moran’s I ≈ +0,28, p ≈ 0,002).

Secara umum, data PM2.5 Provinsi Banten pada studi kasus ini melanggar asumsi independensi yang disyaratkan banyak metode statistik klasik — baik secara temporal maupun spasial — sehingga analisis lanjutan (Pertemuan 3 dan seterusnya: uji kehomogenan dan stasioneritas) perlu memperhitungkan struktur ketergantungan ini agar kesimpulan yang dihasilkan tetap valid.

8 Referensi

  • Barnett, V. (2004). Environmental Statistics: Methods and Applications. Wiley.
  • Rogerson, P. A. (2010). Statistical Methods for Geography. Sage.
  • Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.
  • Tobler, W. R. (1970). A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography, 46, 234–240.
  • Dokumentasi paket R: readxl, dplyr, ggplot2, sf, tmap, spdep (CRAN).

Catatan: seluruh data pada dokumen ini bersifat simulasi/ilustratif untuk keperluan pembelajaran (lihat sheet “Panduan” pada file data). Nilai-nilai pada bagian Interpretasi dihitung dari data simulasi tersebut dengan seed tetap, sehingga akan direproduksi persis apabila dokumen ini di-knit ulang dengan file data yang sama.