1 Pendahuluan

1.1 Tentang Studi Kasus

Dokumen ini menerapkan teknik eksplorasi dan pengujian data lingkungan yang dibahas pada Pertemuan 3 dan 4 (uji kehomogenan dan uji stasioneritas), serta analisis lanjutan berupa ANOVA, uji-t independen, dan forecasting, pada data curah hujan bulanan di wilayah Banten.

Dataset yang digunakan adalah data curah hujan bulanan (mm) dari 4 stasiun pengamatan BMKG di wilayah Banten (Stasiun Klimatologi Tangerang Selatan, Stasiun Meteorologi Serang, Stasiun Meteorologi Curug, dan Stasiun Geofisika Tangerang) sepanjang Januari 2023–Desember 2024 (24 bulan × 4 stasiun = 96 pengamatan).

  • Sumber data: Badan Meteorologi, Klimatologi, dan Geofisika (BMKG) – Data Online. [isi URL halaman data dan tanggal akses, misalnya: https://dataonline.bmkg.go.id, diakses … 2026]
  • Download data: [isi link data (GitHub/Google Drive) di sini]

1.2 Rumusan Pertanyaan Analisis

  1. Bagaimana karakteristik statistik deskriptif dan pola visual curah hujan pada masing-masing stasiun?
  2. Apakah varians curah hujan antarstasiun homogen, dan apakah deret waktunya homogen (tanpa titik perubahan)?
  3. Apakah deret waktu curah hujan stasioner, atau terdapat heterogenitas (tren, ragam yang berubah)?
  4. Apakah rata-rata curah hujan berbeda antarstasiun (ANOVA) dan antartahun (uji-t independen)?
  5. Bagaimana perkiraan curah hujan 12 bulan ke depan (forecasting)?

Kelima pertanyaan ini dijawab bertahap melalui langkah-langkah analisis di bawah ini.

2 Persiapan

2.1 Memuat Package

# Jika paket belum terpasang, jalankan sekali di Console:
# install.packages(c("car", "tseries", "randtests", "trend",
#                    "forecast", "ggplot2", "knitr"))

library(car)        # uji Levene
library(tseries)    # uji ADF dan KPSS
library(randtests)  # Runs test
library(trend)      # uji SNHT
library(forecast)   # forecasting (Holt, ETS, ARIMA)
library(ggplot2)    # visualisasi
library(knitr)      # tabel

2.2 Langkah 1 — Import dan Struktur Data

Nilai - pada data asli diartikan sebagai data tidak tersedia (NA), bukan curah hujan 0 mm. Data ditulis kronologis (2023 lalu 2024).

stasiun_kode <- c("TangSel", "Serang", "Curug", "Geofisika")
stasiun_nama <- c("Klimatologi Tangerang Selatan", "Meteorologi Serang",
                  "Meteorologi Curug", "Geofisika Tangerang")
bulan_id <- c("Januari", "Februari", "Maret", "April", "Mei", "Juni", "Juli",
              "Agustus", "September", "Oktober", "November", "Desember")

data_hujan <- data.frame(
  Tahun = rep(c(2023, 2024), each = 12),
  Bulan = rep(bulan_id, times = 2),
  TangSel = c(110.3, 448.2, 262.4, 76.6, 90.8, 124.9, 53.5, 3.6, 30.4, 13.4, 225.9, 221.5,
              430.5, 212.4, 447, 377.4, 57, 133.6, 202.8, 76.3, 126.6, 25.9, 464.3, 156.4),
  Serang = c(148.7, 310.9, 276.5, 130.1, 69.1, 151.3, 89, NA, 0.5, NA, 138.2, 77.8,
             368, 197.3, 310.5, 163.7, 305.2, 129.5, 58.2, 4.9, 72.2, 92.4, 168.7, 332.8),
  Curug = c(129.9, 246.7, 234.7, 124.9, 78.6, 223.9, 171.4, 78.7, NA, 12.9, 308.6, 156.3,
            322.6, 191, 149.5, 483.6, 102.6, 123.4, 88.7, 234.8, 76.5, 60.1, 613, 218.9),
  Geofisika = c(82.4, 396.4, 278.1, 77.9, 69.7, 155.9, 98.7, 2, 15.5, 24.8, 139.4, 88.3,
                309.7, 160.6, 320.8, 337.3, 169.1, 119.9, 112, 21.8, 104.9, 16.6, 171.4, 294.4)
)

# (Alternatif) membaca langsung dari file Excel asli. Letakkan file di folder yang sama dengan Rmd:
# library(readxl)
# data_hujan <- as.data.frame(read_excel("TUGAS_SL_2.xlsx", na = c("-", "")))
# names(data_hujan) <- c("Tahun", "Bulan", stasiun_kode)
# (lalu urutkan agar 2023 di atas, 2024 di bawah)

data_hujan$Tanggal <- as.Date(sprintf("%d-%02d-01", data_hujan$Tahun,
                                      match(data_hujan$Bulan, bulan_id)))
data_hujan <- data_hujan[order(data_hujan$Tanggal), ]
rownames(data_hujan) <- NULL

# Melihat struktur dan ukuran data
str(data_hujan)
## 'data.frame':    24 obs. of  7 variables:
##  $ Tahun    : num  2023 2023 2023 2023 2023 ...
##  $ Bulan    : chr  "Januari" "Februari" "Maret" "April" ...
##  $ TangSel  : num  110.3 448.2 262.4 76.6 90.8 ...
##  $ Serang   : num  148.7 310.9 276.5 130.1 69.1 ...
##  $ Curug    : num  129.9 246.7 234.7 124.9 78.6 ...
##  $ Geofisika: num  82.4 396.4 278.1 77.9 69.7 ...
##  $ Tanggal  : Date, format: "2023-01-01" "2023-02-01" ...
head(data_hujan[, c("Tahun", "Bulan", stasiun_kode)], 6)
dim(data_hujan)
## [1] 24  7
# Memeriksa data hilang (missing value)
sum(is.na(data_hujan[stasiun_kode]))
## [1] 3
colSums(is.na(data_hujan[stasiun_kode]))
##   TangSel    Serang     Curug Geofisika 
##         0         2         1         0
# Format panjang (long) untuk visualisasi dan uji antarstasiun
long <- data.frame(
  Tanggal = rep(data_hujan$Tanggal, times = 4),
  Tahun   = rep(data_hujan$Tahun,   times = 4),
  Bulan   = rep(data_hujan$Bulan,   times = 4),
  Stasiun = factor(rep(stasiun_nama, each = nrow(data_hujan)), levels = stasiun_nama),
  Nilai   = unlist(data_hujan[stasiun_kode], use.names = FALSE)
)

# Rata-rata wilayah (rata-rata stasiun yang tersedia pada bulan tersebut)
rata_df <- data.frame(Tanggal = data_hujan$Tanggal, Tahun = data_hujan$Tahun,
                      Bulan = data_hujan$Bulan,
                      Rata = rowMeans(data_hujan[stasiun_kode], na.rm = TRUE))
rata_df$Tahun_f <- factor(rata_df$Tahun)

Interpretasi

Data berisi 24 baris dan 7 kolom (termasuk kolom Tanggal yang dibuat untuk pengurutan waktu), yaitu 24 bulan × 4 stasiun. Kolom curah hujan bertipe numerik dan Bulan bertipe karakter. Terdapat 3 nilai hilang, yaitu Serang (Agustus dan Oktober 2023) dan Curug (September 2023). Ketiganya berada di musim kemarau sehingga bisa berarti alat tidak mencatat atau hujan sangat kecil. Karena itu nilai tersebut tidak dianggap 0 mm: pada uji antarstasiun dan ANOVA nilai kosong dikeluarkan otomatis, dan pada uji deret waktu diisi interpolasi linear (dijelaskan di Langkah 7).

3 Statistik Deskriptif

3.1 Langkah 2 — Statistik Deskriptif

# Statistik deskriptif keseluruhan
summary(long$Nilai)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.5    77.9   138.2   171.7   234.8   613.0       3
sd(long$Nilai, na.rm = TRUE)
## [1] 128.6931
IQR(long$Nilai, na.rm = TRUE)
## [1] 156.9
# Statistik deskriptif per stasiun
desk <- do.call(rbind, lapply(split(long$Nilai, long$Stasiun), function(v) {
  data.frame(n = sum(!is.na(v)), Rata2 = mean(v, na.rm = TRUE),
             Median = median(v, na.rm = TRUE),
             SD = sd(v, na.rm = TRUE), Varians = var(v, na.rm = TRUE),
             Min = min(v, na.rm = TRUE), Maks = max(v, na.rm = TRUE))
}))
kable(desk, digits = 1, caption = "Statistik deskriptif curah hujan bulanan (mm) per stasiun")
Statistik deskriptif curah hujan bulanan (mm) per stasiun
n Rata2 Median SD Varians Min Maks
Klimatologi Tangerang Selatan 24 182.2 130.1 149.9 22470.8 3.6 464.3
Meteorologi Serang 22 163.4 143.4 108.9 11863.0 0.5 368.0
Meteorologi Curug 23 192.7 156.3 139.0 19320.3 12.9 613.0
Geofisika Tangerang 24 148.7 116.0 115.3 13291.7 2.0 396.4 
rata_all <- mean(long$Nilai, na.rm = TRUE)
med_all  <- median(long$Nilai, na.rm = TRUE)
st_tinggi <- rownames(desk)[which.max(desk$Rata2)]
st_rendah <- rownames(desk)[which.min(desk$Rata2)]
st_sd_maks <- rownames(desk)[which.max(desk$SD)]

Interpretasi

Secara keseluruhan, rata-rata curah hujan bulanan adalah 171,7 mm dengan median 138,2 mm. Rata-rata yang lebih besar daripada median menandakan distribusi condong ke kanan (right-skewed): beberapa bulan dengan hujan sangat tinggi menarik rata-rata ke atas. Stasiun dengan rata-rata tertinggi adalah Meteorologi Curug dan terendah Geofisika Tangerang, sedangkan simpangan baku terbesar ada pada Klimatologi Tangerang Selatan. Selisih rata-rata antarstasiun tidak besar dibanding simpangan bakunya (variasi antarbulan jauh lebih besar daripada variasi antarstasiun), dan hal ini diuji secara formal pada Langkah 10.

4 Visualisasi Data

4.1 Langkah 3 — Histogram dan Boxplot

# Histogram distribusi curah hujan keseluruhan
ggplot(long, aes(x = Nilai)) +
  geom_histogram(binwidth = 50, fill = "#147D64", color = "white", na.rm = TRUE) +
  labs(title = "Distribusi Curah Hujan Bulanan - Wilayah Banten 2023–2024",
       x = "Curah Hujan (mm)", y = "Frekuensi") +
  theme_minimal()

# Boxplot curah hujan per stasiun
ggplot(long, aes(x = Stasiun, y = Nilai, fill = Stasiun)) +
  geom_boxplot(show.legend = FALSE, na.rm = TRUE) +
  labs(title = "Boxplot Curah Hujan Bulanan Antarstasiun",
       x = "Stasiun", y = "Curah Hujan (mm)") +
  scale_x_discrete(labels = function(x) sub(" ", "\n", x)) +
  theme_minimal()

Interpretasi

Histogram memperlihatkan sebagian besar bulan memiliki curah hujan rendah–menengah dengan ekor memanjang ke kanan menuju bulan-bulan berhujan sangat tinggi (di atas 400 mm), sejalan dengan kecondongan positif pada Langkah 2. Boxplot digunakan sebagai indikasi awal homogenitas varians: jika panjang kotak dan whisker antarstasiun mirip, varians cenderung homogen. Kotak dan whisker keempat stasiun tampak berada pada rentang yang serupa, dan titik-titik di luar whisker (jika ada) menandakan bulan-bulan berhujan ekstrem.

4.2 Langkah 4 — Diagram Deret Waktu

ggplot(long, aes(x = Tanggal, y = Nilai, color = Stasiun)) +
  geom_line(linewidth = 0.9, na.rm = TRUE) + geom_point(size = 1.5, na.rm = TRUE) +
  labs(title = "Curah Hujan Bulanan per Stasiun (2023–2024)",
       x = "Bulan", y = "Curah Hujan (mm)") +
  theme_minimal() + theme(legend.position = "bottom")

# Pola musiman: rata-rata wilayah per bulan, dibandingkan antartahun
rata_df$Bulan_no <- match(rata_df$Bulan, bulan_id)
ggplot(rata_df, aes(x = Bulan_no, y = Rata, color = Tahun_f, group = Tahun_f)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_x_continuous(breaks = 1:12, labels = substr(bulan_id, 1, 3)) +
  labs(title = "Rata-rata Curah Hujan Wilayah per Bulan",
       x = "Bulan", y = "Curah Hujan (mm)", color = "Tahun") +
  theme_minimal()

# Rata-rata musim hujan (Nov–Apr) vs kemarau (Mei–Okt), rata-rata wilayah
rata_df$Musim <- ifelse(rata_df$Bulan_no %in% c(11, 12, 1, 2, 3, 4), "Hujan", "Kemarau")
tapply(rata_df$Rata, rata_df$Musim, mean)
##     Hujan   Kemarau 
## 247.76042  86.12917

Interpretasi

Grafik deret waktu memperlihatkan pola musiman yang konsisten pada keempat stasiun: curah hujan tinggi pada bulan-bulan musim hujan (sekitar November–April) dan menurun tajam pada musim kemarau (sekitar Juli–Oktober), lalu naik kembali. Keempat stasiun bergerak searah, sehingga fluktuasi antarbulan jauh lebih dominan daripada perbedaan antarstasiun. Rata-rata wilayah pada musim hujan sebesar 247,8 mm, sedangkan pada musim kemarau 86,1 mm.

5 Uji Homogenitas Data

5.1 Langkah 5 — Uji Normalitas (Prasyarat)

Bartlett mengasumsikan data normal dan sensitif terhadap pelanggaran normalitas, sedangkan Levene lebih robust. Karena itu normalitas diperiksa dahulu dengan Shapiro-Wilk (H₀: data normal).

sh <- do.call(rbind, lapply(seq_along(stasiun_kode), function(i) {
  s <- shapiro.test(na.omit(data_hujan[[stasiun_kode[i]]]))
  data.frame(Stasiun = stasiun_nama[i], W = unname(s$statistic), p_value = s$p.value,
             Kesimpulan = ifelse(s$p.value < alpha, "Tidak normal", "Normal"))
}))
kable(sh, digits = 4, caption = "Uji normalitas Shapiro-Wilk per stasiun")
Uji normalitas Shapiro-Wilk per stasiun
Stasiun W p_value Kesimpulan
Klimatologi Tangerang Selatan 0.8774 0.0074 Tidak normal
Meteorologi Serang 0.9263 0.1029 Normal
Meteorologi Curug 0.8588 0.0040 Tidak normal
Geofisika Tangerang 0.9077 0.0315 Tidak normal 
n_tn <- sum(sh$p_value < alpha)

Interpretasi

3 dari 4 stasiun menunjukkan data tidak normal pada α = 5%. Hal ini wajar untuk curah hujan bulanan yang menjulur ke kanan. Akibatnya, hasil Levene dijadikan rujukan utama untuk homogenitas varians, sedangkan Bartlett disajikan sebagai pembanding.

5.2 Langkah 6 — Homogenitas Varians Antarstasiun (Bartlett dan Levene)

  • H₀: varians curah hujan semua stasiun sama (homogen)
  • H₁: minimal ada satu stasiun dengan varians berbeda (tidak homogen)
bt <- bartlett.test(Nilai ~ Stasiun, data = long)
lv <- leveneTest(Nilai ~ Stasiun, data = long, center = median)
bt
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Nilai by Stasiun
## Bartlett's K-squared = 2.9642, df = 3, p-value = 0.3972
lv
p_bt <- bt$p.value; K_bt <- unname(bt$statistic); df_bt <- unname(bt$parameter)
p_lv <- lv[["Pr(>F)"]][1]; F_lv <- lv[["F value"]][1]

Interpretasi

Uji Bartlett menghasilkan K² = 2,964 (db = 3) dengan p-value = 0,3972. Uji Levene menghasilkan F = 0,549 dengan p-value = 0,6500. Karena p-value Levene ≥ 0,05, maka gagal tolak H₀: varians curah hujan antarstasiun dapat dianggap homogen (sama). Kesimpulan Bartlett sejalan dengan Levene, dan hal ini sesuai dengan boxplot yang menunjukkan sebaran antarstasiun relatif mirip.

5.3 Langkah 7 — Homogenitas Deret Waktu (SNHT)

Standard Normal Homogeneity Test (SNHT) mendeteksi titik perubahan (break point) rata-rata pada deret waktu.

  • H₀: deret waktu homogen (tidak ada titik perubahan)
  • H₁: terdapat titik perubahan

Nilai kosong pada tiap stasiun diisi dengan interpolasi linear hanya untuk uji deret waktu (Langkah 7 sampai 9 dan forecasting memakai rata-rata wilayah tanpa imputasi), agar deret tidak terputus.

set.seed(123)
mk_ts  <- function(x) ts(x, start = c(2023, 1), frequency = 12)
interp <- function(x) {
  idx <- seq_along(x); ok <- !is.na(x)
  approx(idx[ok], x[ok], xout = idx, rule = 2)$y
}

seri_list <- c(lapply(data_hujan[stasiun_kode], function(x) mk_ts(interp(x))),
               list(RataRata = mk_ts(rata_df$Rata)))
nama_seri <- c(stasiun_nama, "Rata-rata wilayah")

snht_tab <- do.call(rbind, lapply(seq_along(seri_list), function(i) {
  s <- trend::snh.test(seri_list[[i]])
  k <- as.integer(unname(s$estimate))
  data.frame(Seri = nama_seri[i], Statistik_T = unname(s$statistic),
             p_value = s$p.value, K = k,
             Titik_ubah = paste(data_hujan$Bulan[k], data_hujan$Tahun[k]),
             Kesimpulan = ifelse(s$p.value < alpha, "Tidak homogen", "Homogen"))
}))
kable(snht_tab, digits = 4, caption = "Hasil SNHT per deret waktu")
Hasil SNHT per deret waktu
Seri Statistik_T p_value K Titik_ubah Kesimpulan
Klimatologi Tangerang Selatan 2.8150 0.6437 10 Oktober 2023 Homogen
Meteorologi Serang 2.8219 0.6396 23 November 2024 Homogen
Meteorologi Curug 5.9250 0.1416 22 Oktober 2024 Homogen
Geofisika Tangerang 2.7712 0.6514 3 Maret 2023 Homogen
Rata-rata wilayah 3.2045 0.5502 22 Oktober 2024 Homogen 
n_snht <- sum(snht_tab$p_value < alpha)

Interpretasi

Berdasarkan SNHT, 0 dari 5 deret (seluruhnya) tidak menunjukkan titik perubahan yang signifikan pada α = 5%. Artinya deret waktu 2023–2024 dapat dianggap homogen sepanjang periode. Dengan hanya 24 titik, kekuatan uji juga terbatas.

6 Uji Heterogenitas: Stasioneritas Data

6.1 Langkah 8 — Uji Stasioneritas (ADF, KPSS, Runs Test, ARCH-LM)

Uji H₀ Kriteria (α = 5%)
ADF ada unit root (tidak stasioner) p < 0,05 → stasioner
KPSS (level) data stasioner p < 0,05 → tidak stasioner
Runs Test data acak p < 0,05 → tidak acak (ada pola)
ARCH-LM tidak ada efek ARCH (ragam konstan) p < 0,05 → ragam tidak konstan

Uji ARCH-LM dihitung manual dengan meregresikan kuadrat deviasi dari rata-rata pada 2 lag sebelumnya. Statistiknya (n − q)·R² menyebar χ² dengan db = q.

par(mfrow = c(1, 2))
plot(seri_list$RataRata, type = "b", col = "darkgreen", pch = 19,
     main = "Rata-rata Curah Hujan Wilayah", ylab = "mm", xlab = "Waktu")
acf(seri_list$RataRata, main = "ACF Rata-rata Wilayah", lag.max = 12)

par(mfrow = c(1, 1))
arch_lm <- function(x, q = 2) {
  e2 <- (x - mean(x))^2
  n  <- length(e2)
  d  <- embed(e2, q + 1)
  R2 <- summary(lm(d[, 1] ~ d[, -1]))$r.squared
  pchisq((n - q) * R2, df = q, lower.tail = FALSE)
}

uji_stasioner <- function(x) {
  x <- as.numeric(x)
  adf <- suppressWarnings(tseries::adf.test(x))
  kp  <- suppressWarnings(tseries::kpss.test(x, null = "Level"))
  rn  <- randtests::runs.test(x)
  data.frame(ADF_p = adf$p.value, KPSS_p = kp$p.value,
             Runs_p = rn$p.value, ARCH_p = arch_lm(x, 2))
}

st <- do.call(rbind, lapply(seri_list, uji_stasioner))
st <- data.frame(Seri = nama_seri, st, row.names = NULL)
st$ADF   <- ifelse(st$ADF_p  < alpha, "Stasioner", "Tidak stasioner")
st$KPSS  <- ifelse(st$KPSS_p < alpha, "Tidak stasioner", "Stasioner")
st$Runs  <- ifelse(st$Runs_p < alpha, "Tidak acak", "Acak")
st$ARCH  <- ifelse(st$ARCH_p < alpha, "Ragam tidak konstan", "Ragam konstan")

kable(st[, c("Seri", "ADF_p", "ADF", "KPSS_p", "KPSS", "Runs_p", "Runs", "ARCH_p", "ARCH")],
      digits = 4, caption = "Ringkasan uji stasioneritas (nilai p dan kesimpulan)")
Ringkasan uji stasioneritas (nilai p dan kesimpulan)
Seri ADF_p ADF KPSS_p KPSS Runs_p Runs ARCH_p ARCH
Klimatologi Tangerang Selatan 0.4362 Tidak stasioner 0.1 Stasioner 0.0369 Tidak acak 0.7671 Ragam konstan
Meteorologi Serang 0.3292 Tidak stasioner 0.1 Stasioner 0.0950 Acak 0.2845 Ragam konstan
Meteorologi Curug 0.3959 Tidak stasioner 0.1 Stasioner 0.6764 Acak 0.8376 Ragam konstan
Geofisika Tangerang 0.4324 Tidak stasioner 0.1 Stasioner 0.2105 Acak 0.6714 Ragam konstan
Rata-rata wilayah 0.3871 Tidak stasioner 0.1 Stasioner 0.0369 Tidak acak 0.4491 Ragam konstan 
n_adf  <- sum(st$ADF  == "Stasioner")
n_kpss <- sum(st$KPSS == "Stasioner")
n_runs <- sum(st$Runs == "Acak")
n_arch <- sum(st$ARCH == "Ragam konstan")

Catatan: p-value KPSS pada paket tseries dibatasi pada rentang 0,01–0,10, sehingga nilai 0,01 berarti ≤ 0,01 dan 0,10 berarti ≥ 0,10.

Interpretasi

  • ADF: 0 dari 5 deret dinyatakan stasioner (tolak H₀ unit root).
  • KPSS: 5 dari 5 deret dinyatakan stasioner (gagal tolak H₀).
  • Runs Test: 3 dari 5 deret dinyatakan acak. Uji ini menguji keacakan, bukan stasioneritas secara langsung: deret yang stasioner tetap bisa tidak acak jika memiliki autokorelasi atau pola musiman.
  • ARCH-LM: 5 dari 5 deret memiliki ragam yang konstan (tidak ada efek heteroskedastisitas).

Hasil ADF dan KPSS belum sepenuhnya seragam antarderet; bila ada deret yang tidak stasioner terhadap rataan, penanganan yang sesuai adalah differencing. Dengan hanya 24 pengamatan (dua siklus musiman), uji ADF berdaya rendah dan pola musiman belum dapat dipisahkan dari fluktuasi acak secara meyakinkan.

6.2 Langkah 9 — Homogenitas Varians Antartahun (2023 vs 2024)

lv_th <- leveneTest(Rata ~ Tahun_f, data = rata_df)
lv_th
p_lv_th <- lv_th[["Pr(>F)"]][1]

Interpretasi

Uji Levene pada rata-rata wilayah menghasilkan p-value = 0,3844, sehingga varians curah hujan antara 2023 dan 2024 dapat dianggap homogen.

7 Analisis Lanjutan

7.1 Langkah 10 — ANOVA Satu Arah: Perbedaan Rata-rata Antarstasiun

  • H₀: μ₁ = μ₂ = μ₃ = μ₄ (rata-rata curah hujan semua stasiun sama)
  • H₁: minimal ada satu stasiun dengan rata-rata berbeda
fit_aov <- aov(Nilai ~ Stasiun, data = long)
summary(fit_aov)
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Stasiun      3   26988    8996   0.535   0.66
## Residuals   89 1496709   16817               
## 3 observations deleted due to missingness
tab_aov <- summary(fit_aov)[[1]]
F_aov <- tab_aov[["F value"]][1]; p_aov <- tab_aov[["Pr(>F)"]][1]
db1 <- tab_aov[["Df"]][1]; db2 <- tab_aov[["Df"]][2]

sh_res <- shapiro.test(residuals(fit_aov))
kw <- kruskal.test(Nilai ~ Stasiun, data = long)
sh_res
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit_aov)
## W = 0.92293, p-value = 3.822e-05
kw
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Nilai by Stasiun
## Kruskal-Wallis chi-squared = 1.2769, df = 3, p-value = 0.7346
aggregate(Nilai ~ Stasiun, data = long, FUN = mean)
TukeyHSD(fit_aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Nilai ~ Stasiun, data = long)
## 
## $Stasiun
##                                                        diff        lwr
## Meteorologi Serang-Klimatologi Tangerang Selatan  -18.72235 -118.94025
## Meteorologi Curug-Klimatologi Tangerang Selatan    10.51105  -88.56364
## Geofisika Tangerang-Klimatologi Tangerang Selatan -33.50417 -131.51921
## Meteorologi Curug-Meteorologi Serang               29.23340  -72.02110
## Geofisika Tangerang-Meteorologi Serang            -14.78182 -114.99972
## Geofisika Tangerang-Meteorologi Curug             -44.01522 -143.08991
##                                                         upr     p adj
## Meteorologi Serang-Klimatologi Tangerang Selatan   81.49555 0.9613292
## Meteorologi Curug-Klimatologi Tangerang Selatan   109.58574 0.9924652
## Geofisika Tangerang-Klimatologi Tangerang Selatan  64.51087 0.8074584
## Meteorologi Curug-Meteorologi Serang              130.48790 0.8739009
## Geofisika Tangerang-Meteorologi Serang             85.43608 0.9803038
## Geofisika Tangerang-Meteorologi Curug              55.05947 0.6514967

Interpretasi

Uji asumsi. Homogenitas varians terpenuhi jika p-value Levene ≥ 0,05 (p = 0,6500). Normalitas residual diuji dengan Shapiro-Wilk (p = < 0,001), sehingga residual tidak normal. Oleh karena itu Kruskal-Wallis sebagai alternatif nonparametrik ikut disajikan.

Hasil. ANOVA menghasilkan F(3; 89) = 0,535 dengan p-value = 0,6595. Karena p-value ≥ 0,05, gagal tolak H₀: tidak terdapat perbedaan rata-rata curah hujan yang signifikan antarstasiun. Kruskal-Wallis memberikan p-value = 0,7346, yang menguatkan hasil ANOVA. Uji lanjut Tukey HSD ditampilkan sebagai pelengkap: bila ANOVA tidak signifikan, seluruh pasangan stasiun seharusnya juga tidak berbeda nyata (p adj ≥ 0,05).

Catatan keterbatasan. Keempat stasiun diamati pada bulan yang sama sehingga datanya berpasangan menurut waktu, sedangkan ANOVA satu arah mengasumsikan kelompok saling bebas. Hasil ini sebaiknya dibaca sebagai analisis eksploratif.

7.2 Langkah 11 — Uji-t Independen: 2023 vs 2024

Uji-t membandingkan rata-rata curah hujan bulanan wilayah (rata-rata 4 stasiun) antara tahun 2023 dan 2024 (n = 12 bulan per tahun).

  • H₀: μ₂₀₂₃ = μ₂₀₂₄
  • H₁: μ₂₀₂₃ ≠ μ₂₀₂₄
m_th <- tapply(rata_df$Rata, rata_df$Tahun_f, mean)
sh_th <- tapply(rata_df$Rata, rata_df$Tahun_f, function(v) shapiro.test(v)$p.value)
tt <- t.test(Rata ~ Tahun_f, data = rata_df)   # Welch (tidak mengasumsikan varians sama)
wt <- wilcox.test(Rata ~ Tahun_f, data = rata_df, exact = FALSE)
m_th
##     2023     2024 
## 131.4563 202.4333
sh_th
##      2023      2024 
## 0.3079140 0.1442092
tt
## 
##  Welch Two Sample t-test
## 
## data:  Rata by Tahun_f
## t = -1.6071, df = 21.696, p-value = 0.1225
## alternative hypothesis: true difference in means between group 2023 and group 2024 is not equal to 0
## 95 percent confidence interval:
##  -162.64438   20.69021
## sample estimates:
## mean in group 2023 mean in group 2024 
##           131.4563           202.4333
wt
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Rata by Tahun_f
## W = 47, p-value = 0.1572
## alternative hypothesis: true location shift is not equal to 0

Interpretasi

Rata-rata curah hujan bulanan wilayah pada 2023 sebesar 131,5 mm dan pada 2024 sebesar 202,4 mm. Uji-t Welch menghasilkan t = -1,607, db = 21,7, p-value = 0,1225. Karena p-value ≥ 0,05, gagal tolak H₀: tidak ada perbedaan rata-rata yang signifikan antara 2023 dan 2024. Uji Wilcoxon (nonparametrik) memberikan p-value = 0,1572, yang sejalan dengan uji-t. Dengan hanya 12 bulan per tahun, daya uji terbatas sehingga selisih rata-rata yang tampak besar belum tentu signifikan.

8 Forecasting

Peramalan dilakukan pada rata-rata curah hujan wilayah (rata-rata 4 stasiun, tanpa imputasi) untuk 12 bulan ke depan (Januari–Desember 2025).

8.1 Langkah 12 — Pemilihan Model dengan Data Uji

Data 24 bulan dibagi menjadi data latih (Jan 2023 – Jun 2024, 18 bulan) dan data uji (Jul – Des 2024, 6 bulan). Model yang dibandingkan adalah Mean, Seasonal naive, Holt (double exponential smoothing, seperti pada materi), ETS, dan ARIMA (keduanya non-musiman karena data belum cukup untuk estimasi musiman yang andal).

rata_ts <- mk_ts(rata_df$Rata)
train <- window(rata_ts, end = c(2024, 6))
test  <- window(rata_ts, start = c(2024, 7))

mod_fun <- list(
  "Mean"                = function(y, h) meanf(y, h = h),
  "Seasonal naive"      = function(y, h) snaive(y, h = h),
  "Holt"                = function(y, h) holt(y, h = h),
  "ETS (non-musiman)"   = function(y, h) forecast(ets(y, model = "ZZN"), h = h),
  "ARIMA (non-musiman)" = function(y, h) forecast(auto.arima(y, seasonal = FALSE), h = h)
)

hasil <- do.call(rbind, lapply(names(mod_fun), function(nm) {
  f <- mod_fun[[nm]](train, h = length(test))
  e <- as.numeric(test) - as.numeric(f$mean)
  data.frame(Model = nm, RMSE = sqrt(mean(e^2)), MAE = mean(abs(e)))
}))
hasil <- hasil[order(hasil$RMSE), ]
kable(hasil, digits = 2, row.names = FALSE,
      caption = "Akurasi peramalan pada data uji (Jul–Des 2024)")
Akurasi peramalan pada data uji (Jul–Des 2024)
Model RMSE MAE
Seasonal naive 88.23 74.32
ETS (non-musiman) 108.85 100.18
ARIMA (non-musiman) 108.85 100.18
Mean 108.85 100.18
Holt 110.70 103.24 
best <- hasil$Model[1]

Interpretasi

Model dengan RMSE terkecil pada data uji adalah Seasonal naive (RMSE = 88,2 mm). Model sederhana justru terbaik, sehingga model yang lebih kompleks tidak memberikan keunggulan berarti pada data sependek ini.

8.2 Langkah 13 — Peramalan 12 Bulan ke Depan

Model terpilih dilatih ulang pada seluruh 24 bulan data.

fc_final <- mod_fun[[best]](rata_ts, 12)

plot(fc_final, main = paste("Forecast Curah Hujan Wilayah –", best),
     ylab = "Curah Hujan (mm)", xlab = "Tahun", col = "darkred")

tab_fc <- data.frame(
  Bulan = paste(bulan_id, 2025),
  Prediksi = as.numeric(fc_final$mean),
  Bawah_95 = pmax(as.numeric(fc_final$lower[, 2]), 0),
  Atas_95  = as.numeric(fc_final$upper[, 2])
)
kable(tab_fc, digits = 1, caption = "Prediksi curah hujan rata-rata wilayah (mm) beserta selang 95%")
Prediksi curah hujan rata-rata wilayah (mm) beserta selang 95%
Bulan Prediksi Bawah_95 Atas_95
Januari 2025 357.7 106.8 608.6
Februari 2025 190.3 0.0 441.2
Maret 2025 307.0 56.1 557.8
April 2025 340.5 89.6 591.4
Mei 2025 158.5 0.0 409.3
Juni 2025 126.6 0.0 377.5
Juli 2025 115.4 0.0 366.3
Agustus 2025 84.4 0.0 335.3
September 2025 95.0 0.0 345.9
Oktober 2025 48.8 0.0 299.6
November 2025 354.4 103.5 605.2
Desember 2025 250.6 0.0 501.5

Interpretasi

Nilai prediksi disertai selang kepercayaan 95% (batas bawah dipotong pada 0 mm karena curah hujan tidak bisa negatif). Selang yang lebar menandakan ketidakpastian tinggi. Dengan hanya dua tahun data yang bermusim kuat, model non-musiman cenderung menghasilkan prediksi yang datar dan belum menangkap siklus musim hujan–kemarau. Karena itu hasil ini bersifat indikatif. Untuk peramalan yang andal diperlukan deret yang lebih panjang (idealnya ≥ 10 tahun) agar model musiman (misalnya SARIMA atau Holt-Winters) dapat diestimasi.

9 Kesimpulan

Menjawab kelima pertanyaan analisis di awal:

  1. Deskriptif dan visualisasi: rata-rata curah hujan bulanan tertinggi ada di Meteorologi Curug dan terendah di Geofisika Tangerang. Distribusi condong ke kanan, dan terlihat pola musiman yang jelas (tinggi pada musim hujan, rendah pada Juli–Oktober) yang searah di keempat stasiun.
  2. Homogenitas: uji Levene (p = 0,6500) dan Bartlett (p = 0,3972) menyimpulkan varians antarstasiun homogen. SNHT menunjukkan 0 dari 5 deret dengan titik perubahan signifikan .
  3. Heterogenitas (stasioneritas): ADF menyatakan 0 dari 5 deret stasioner, KPSS 5 dari 5 deret, dan ARCH-LM menyatakan 5 dari 5 deret berragam konstan.
  4. Perbandingan rata-rata: ANOVA menunjukkan rata-rata antarstasiun tidak berbeda nyata (p = 0,6595), dan uji-t menunjukkan rata-rata 2023 dan 2024 tidak berbeda nyata (p = 0,1225).
  5. Forecasting: model terbaik pada data uji adalah Seasonal naive. Hasil peramalan bersifat indikatif karena keterbatasan panjang data.

Secara umum, data curah hujan ini bermusim kuat dan tidak normal, sehingga uji yang bersifat robust dan nonparametrik dilengkapi bersama uji parametrik.

Keterbatasan: hanya 24 bulan data, ada 3 data kosong (Serang Agustus dan Oktober 2023; Curug September 2023), distribusi tidak normal, dan pola musiman kuat. Kesimpulan uji stasioneritas dan peramalan perlu ditafsirkan hati-hati.

10 Referensi