Pendahuluan

Deskripsi Data

Data yang digunakan pada tugas ini adalah data curah hujan dekadal (10-harian) tingkat provinsi di Indonesia periode Januari 2022 - Agustus 2026, beserta ringkasan tahunannya, yang dilengkapi dengan data luas kebakaran hutan dan lahan (karhutla) per provinsi tahun 2022-2026.

Variabel penting:

  • rfh : curah hujan 10-harian (mm)
  • rfh_avg : rata-rata jangka panjang curah hujan 10-harian (mm)
  • rfq : anomali curah hujan (%), yaitu rasio curah hujan aktual terhadap rata-rata jangka panjang

Tujuan

  1. Menguji homogenitas (kesamaan ragam) curah hujan antar tahun, antara Jawa dan Luar Jawa, antar 6 kelompok pulau, antar 34 provinsi, serta pada provinsi Kalimantan Barat (Kalbar).
  2. Menguji stasioneritas data curah hujan dekadal (deret waktu) Kalbar dengan uji ADF dan KPSS, serta memeriksa pengaruh differencing.
  3. Membandingkan rata-rata curah hujan antar tahun, antar pulau, antar provinsi, dan antara Jawa vs Luar Jawa menggunakan ANOVA (dengan uji lanjut Tukey), Kruskal-Wallis, dan uji-t independen.
  4. Meramalkan curah hujan dekadal Kalbar satu tahun ke depan (36 dekad) dengan model ARIMA.
  5. Menginterpretasikan hasil seluruh uji dan menarik kesimpulan.

Paket yang digunakan

library(readxl)
library(dplyr)
library(ggplot2)
library(car)
library(tseries)
library(forecast)

1. Import dan Persiapan Data

path_excel <- "C:/Users/user/OneDrive - untirta.ac.id/Documents/Semester 5/Datacoba1 (1) FIX.xlsx"

dp <- read_excel(path_excel, sheet="Data_Project") %>%
  mutate(rfh=as.numeric(rfh), tahun=as.factor(tahun)) %>% filter(!is.na(rfh))

1.1 Jawa-Luar jawa

prov_jawa <- c("DKI Jakarta", "Jawa Barat", "Jawa Tengah", "Jawa Timur", "Banten", "Yogyakarta", "DI Yogyakarta")

dp <- dp %>% mutate(
  region = ifelse(provinsi %in% prov_jawa, "Jawa", "Luar Jawa"),
  pulau_6 = case_when(
    provinsi %in% c("Aceh","Sumatera Utara","Sumatera Barat","Riau","Jambi","Sumatera Selatan","Bengkulu","Lampung","Bangka Belitung","Kepulauan Riau") ~ "Sumatera",
    provinsi %in% c("DKI Jakarta","Jawa Barat","Jawa Tengah","Yogyakarta","Jawa Timur","Banten","DI Yogyakarta") ~ "Jawa",
    provinsi %in% c("Bali","Nusa Tenggara Barat","Nusa Tenggara Timur") ~ "Bali-Nusa",
    provinsi %in% c("Kalimantan Barat","Kalimantan Tengah","Kalimantan Selatan","Kalimantan Timur","Kalimantan Utara") ~ "Kalimantan",
    provinsi %in% c("Sulawesi Utara","Sulawesi Tengah","Sulawesi Selatan","Sulawesi Tenggara","Gorontalo","Sulawesi Barat") ~ "Sulawesi",
    TRUE ~ "Maluku-Papua"
  ),
  is_kalbar = ifelse(provinsi == "Kalimantan Barat", "Kalbar", "Non-Kalbar")
)

# Print cek Jawa
dp %>% filter(region=="Jawa") %>% distinct(provinsi) %>% arrange(provinsi)

Pembagian Region membagi 34 provinsi yaitu Jawa: 6 provinsi (Banten, DKI Jakarta, Jawa Barat, Jawa Tengah, Jawa Timur, Yogyakarta). Luar Jawa: 28 provinsi sisanya.

1.2 Kalbar dekadal

## 1.2 Kalbar dekadal
ch <- read_excel(path_excel, sheet="Curah Hujan 2022-2026")
kalbar <- ch %>% filter(adm_level==1, PCODE=="ID61") %>% arrange(date) %>% select(date, rfh, rfh_avg, rfq)
y <- ts(kalbar$rfh, start=c(2022,1), frequency=36)

head(kalbar)
print(y)
## Time Series:
## Start = c(2022, 1) 
## End = c(2026, 24) 
## Frequency = 36 
##   [1] 118.23538 156.28317  88.68308 192.53218  84.74052  39.84489  95.74806
##   [8] 122.11382  90.83232 101.59652 104.84531  77.25404 113.86586 139.97128
##  [15]  66.10187 102.96185 146.11444 109.35066  52.91029  55.41354 133.53406
##  [22]  91.84469  69.79145 216.71913 117.62859 142.60553 111.55649 179.46678
##  [29] 138.98532 138.84384 145.84846 114.91888 102.43178 152.19011 127.55376
##  [36]  86.80885  78.67806  76.43052 181.97778  90.12052  98.03040 137.52127
##  [43] 130.31796 167.26263 142.70530 111.88850  59.03626 102.41585  93.48606
##  [50]  53.04255  66.15887  85.56906  94.06623  44.84092 116.85894  65.26074
##  [57]  25.27185  73.96227  51.37581  60.94173  37.77678  81.98344  49.88702
##  [64]  67.33997  93.81094 126.85558 106.00545 111.32153 139.22113 148.38838
##  [71] 108.09537 136.73590 141.57410 173.59862  95.68162  88.17837 170.20792
##  [78]  87.41899 169.89163  72.10396 152.54265  97.27457 118.17586 116.95661
##  [85] 170.60114 139.64200 123.72312  81.93691 100.13017  76.82414 124.51247
##  [92]  19.79459  17.67156  84.39069  85.68015 151.17564  85.01237  13.86271
##  [99] 125.41627 159.94885 141.03186  61.87403 112.38357 119.32739 133.52043
## [106] 103.21484 145.10773  92.26514 148.93188 159.83882 176.72961  66.58835
## [113]  95.41396 115.30853 151.47935 160.85768 113.29365 104.86397 106.58583
## [120]  89.23224  94.14022  92.99916 101.31335 116.35045  96.73632  74.34794
## [127]  55.04737  58.90547  24.51100 123.42423 145.65270  85.60805 134.45903
## [134] 137.62859  41.49885  63.58520 147.54768  92.50828  95.91972 113.58541
## [141] 130.28757 124.77426 138.42989 182.70992 137.71265  35.57912  61.96898
## [148]  88.19744 153.37749  87.88870 120.13037  77.22385  81.67575 122.72794
## [155]  77.65437  98.34396  93.99183 120.43932  96.73758  61.91365  85.77678
## [162]  60.61224  19.53511  22.27353  35.61916  12.59610  16.91574  20.70509

Data tahunan mencakup 170 observasi (provinsi x tahun, 2022–2026), sedangkan data dekadal Kalbar mencakup 168 periode 10-harian dari Januari 2022 hingga Agustus 2026.

1.3 Summary 34 provinsi

dp_34 <- dp %>% group_by(provinsi) %>% summarise(mean_rfh=mean(rfh), sd_rfh=sd(rfh), n=n()) %>% arrange(mean_rfh)

dp_34

Tabel tersebut merangkum mean_rfh = rata-rata curah hujan, sd_rfh = standar deviasi / keragaman per tahun, dan n = 5 = jumlah tahun pengamatan 2022-2026 untuk 34 provinsi.

Secara umum, terlihat HETEROGENITAS spasial yang sangat jelas. Rentang mean hanya untuk 10 provinsi terendah saja sudah 52.64 mm sampai 70.87 mm.

  1. Nusa Tenggara Timur - 52.64 ± 7.81 mm
    Rata-rata terendah se-Indonesia. Nilai SD paling kecil di antara 10 ini. Artinya: NTT tidak hanya paling kering, tapi juga paling konsisten kering setiap tahunnya. Ini karakteristik iklim semi-arid. Secara statistik ini adalah outlier.

  2. Gorontalo - 55.32 ± 11.44 mm dan NTB - 57.44 ± 8.15 mm
    Klaster kering kedua. NTB masih satu region iklim dengan NTT Bali-Nusa Tenggara. Gorontalo unik karena di Sulawesi tapi jauh lebih kering dari Sulawesi lainnya, indikasi ada anomali lokal.

  3. DKI Jakarta - 57.69 ± 9.24 mm
    Satu-satunya provinsi di Jawa bagian barat yang masuk 5 besar terkering. Ini tidak sesuai pola monsunal Jawa yang seharusnya basah. Bisa di pengaruhi efek urban atau lainnya

  4. Bali 63.68 ± 11.76, Jawa Timur 64.72 ± 12.34, Sulawesi Tengah 65.60 ± 12.23, Yogyakarta 65.88 ± 14.41
    Ini adalah kelompok transisi Jawa Timur - Bali. Mean naik, tapi SD juga naik di atas 11. Terutama Yogyakarta SD 14.41 artinya curah hujannya tidak stabil antar tahun, kadang basah kadang sangat kering.

  5. Sulawesi Barat 70.84 ± 16.02 & Sulawesi Tenggara 70.87 ± 15.85
    Mean sudah mendekati 70 mm, tapi punya SD TERBESAR di top 10 ini (>15.8) Koefisien Variasi (CV = sd/mean) nya >22%. Artinya: paling heterogen secara temporal. Walaupun rata-ratanya rendah, fluktuasi tahunannya sangat tinggi.

2. Eksplorasi Data (Visualisasi)

2.1 Tahunan Antar Tahun

ggplot(dp, aes(x = tahun, y = rfh, fill = tahun)) +
  geom_boxplot(show.legend = FALSE) +
  labs(title = "Boxplot Curah Hujan Rata-rata Tahunan (rfh) Seluruh Provinsi",
       subtitle = "Tiap titik = satu provinsi, dikelompokkan per tahun",
       x = "Tahun", y = "Curah Hujan 10-harian rata-rata (mm)") +
  theme_minimal()

Boxplot ini membandingkan sebaran rata-rata curah hujan 10-harian 34 provinsi di setiap tahun. Dari gambar terlihat:

Median berbeda antar tahun. 2022 ≈ 90 mm, 2023 ≈ 67 mm, 2024 ≈ 80 mm, 2025 ≈ 95 mm (tertinggi), dan 2026 ≈ 65 mm (terendah). Untuk 2023 dan 2026 tergolong kering. Kotak dan mediannya paling rendah. Untuk 2023, ini kemungkinan terkait kondisi kering (misalnya pengaruh El Niño). Untuk 2026, perlu hati-hati karena data baru sampai Agustus, sehingga rata-ratanya belum mencakup Sep-Des yang di banyak wilayah adalah musim hujan. 2022 dan 2025 tergolong basah. Median dan kotaknya paling tinggi, dengan nilai maksimum 2025 mendekati 120 mm. Lebar kotak (IQR) dan whisker relatif mirip, sekitar 15-30 mm di semua tahun. Ini mendukung dugaan awal ragam yang homogen, dan sejalan dengan uji Bartlett (p = 0.105) dan Levene (p = 0.114). Ada satu pencilan (outlier) di 2026 sekitar 94 mm, yaitu provinsi yang jauh lebih basah dibanding provinsi lain pada tahun itu.

2.2 Curah Hujan Tahunan: Jawa vs Luar Jawa

ggplot(dp, aes(x=region, y=rfh, fill=region)) +
  geom_boxplot(show.legend=FALSE) +
  labs(title="2.2 Jawa vs Luar Jawa", y="rfh (mm)") + theme_minimal()

Interpretasi: panjang kotak (IQR) dan whisker pada kedua boxplot terlihat relatif mirip antar kelompok, mengindikasikan dugaan awal ragam yang cenderung homogen -akan dikonfirmasi dengan uji statistik formal pada bagian berikutnya.

2.3 Kalbar Dekadal

ggplot(kalbar, aes(x = date, y = rfh)) +
  geom_line(color = "steelblue") +
  geom_point(size = 0.8, color = "steelblue") +
  labs(title = "Curah Hujan Dekadal Provinsi Kalbar (2022-2026)",
       x = "Tanggal", y = "Curah Hujan 10-harian (mm)") +
  theme_minimal()

Curah hujan dekadal Kalbar berfluktuasi di sekitar rata-rata ≈ 100 mm tanpa tren jangka panjang yang jelas, dengan penurunan berulang pada pertengahan tahun yang menunjukkan pola musiman lemah. Pola ini mendukung dugaan bahwa data cenderung stasioner terhadap rata-rata.

2.4 Per 6 Pulau Besar (Sumatera, Jawa, Bali-Nusa, Kalimantan, Sulawesi, Maluku-Papua)

ggplot(dp, aes(x=reorder(pulau_6, rfh, median), y=rfh, fill=pulau_6)) +
  geom_boxplot(show.legend=FALSE) + 
  geom_jitter(width=0.2, alpha=0.3) +
  coord_flip() +
  labs(title="2.4 Boxplot per 6 Pulau - Bali-Nusa vs Kalimantan",
       x="Pulau", y="rfh (mm)") + 
  theme_minimal()

dp %>% group_by(pulau_6) %>% summarise(mean=mean(rfh)) %>% arrange(mean)

Boxplot menunjukkan perbedaan median curah hujan yang jelas antar pulau, dengan Bali-Nusa paling kering dan Kalimantan paling basah. Lebar sebaran juga berbeda antar pulau (Maluku-Papua dan Sulawesi paling lebar), sehingga ragam tidak homogen

2.5 Per 34 Provinsi ( Top 5 Terkering vs Terbasah)

# Barplot 34 provinsi 
ggplot(dp_34, aes(x=reorder(provinsi, mean_rfh), y=mean_rfh)) +
  geom_col(fill="steelblue") +
  coord_flip() +
  labs(title="2.5 Rata-rata Curah Hujan per 34 Provinsi (Full)",
       x="", y="mean rfh (mm)") + theme_minimal()

# Top 5 terkering vs terbasah
dp_top <- bind_rows(
  dp_34 %>% head(5) %>% mutate(kategori="Terkering"),
  dp_34 %>% tail(5) %>% mutate(kategori="Terbasah")
)
ggplot(dp_top, aes(x=reorder(provinsi, mean_rfh), y=mean_rfh, fill=kategori)) +
  geom_col() + geom_text(aes(label=round(mean_rfh,1)), hjust=-0.1) +
  coord_flip() + ylim(0,110) +
  labs(title="2.5b Top 5 Terkering vs Terbasah dari 34 Provinsi",
       x="", y="mean rfh (mm)") + theme_minimal()

3.Uji Homogenitas Data

Uji homogenitas menguji apakah ragam (varians) curah hujan sama antar kelompok. Hipotesis umum:

  • \(H_0\) : ragam semua kelompok sama (homogen)
  • \(H_1\) : minimal ada satu kelompok dengan ragam berbeda (tidak homogen)

3.1 Antar Tahun

bartlett.test(rfh ~ tahun, data=dp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by tahun
## Bartlett's K-squared = 7.6639, df = 4, p-value = 0.1047
leveneTest(rfh ~ tahun, data=dp)

3.2 Jawa vs Luar Jawa

bartlett.test(rfh ~ region, data=dp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by region
## Bartlett's K-squared = 0.321, df = 1, p-value = 0.571
leveneTest(rfh ~ region, data=dp)

3.3 Antar 6 Pulau

bartlett.test(rfh ~ pulau_6, data=dp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by pulau_6
## Bartlett's K-squared = 14.5, df = 5, p-value = 0.01273
leveneTest(rfh ~ pulau_6, data=dp)
oneway.test(rfh ~ pulau_6, data=dp, var.equal=FALSE) # Welch jika tidak homogen
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  rfh and pulau_6
## F = 18.545, num df = 5.000, denom df = 62.548, p-value = 2.97e-11

3.4 Antar 34 Provinsi

bartlett.test(rfh ~ provinsi, data=dp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by provinsi
## Bartlett's K-squared = 23.12, df = 33, p-value = 0.8997
leveneTest(rfh ~ provinsi, data=dp)

3.5 Homogenitas KHUSUS Kalbar

kalbar_5th <- dp %>% filter(provinsi=="Kalimantan Barat")


#DATA DEKADAL 168 BARIS
kalbar_tahunan <- kalbar %>% 
  mutate(tahun = as.factor(format(date, "%Y")))

# Apakah curah hujan Kalbar stabil antar tahun? (pakai 168 dekad)
bartlett.test(rfh ~ tahun, data=kalbar_tahunan)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by tahun
## Bartlett's K-squared = 0.64915, df = 4, p-value = 0.9574
leveneTest(rfh ~ tahun, data=kalbar_tahunan)
# Apakah Kalbar vs Non-Kalbar homogen? (ini pakai dp yang 170 baris, aman)
bartlett.test(rfh ~ is_kalbar, data=dp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  rfh by is_kalbar
## Bartlett's K-squared = 0.075558, df = 1, p-value = 0.7834
leveneTest(rfh ~ is_kalbar, data=dp)
library(dplyr)
library(car)
library(knitr)

uji_homogen <- function(label, formula, data, alpha = 0.05) {
  b <- bartlett.test(formula, data = data)$p.value
  l <- leveneTest(formula, data = data)$`Pr(>F)`[1]
  tibble(
    `Sub-bagian`  = label,
    `Bartlett p`  = round(b, 3),
    `Levene p`    = round(l, 3),
    Keputusan     = ifelse(b > alpha & l > alpha,
                           "H0 diterima, homogen",
                           ifelse(b <= alpha & l <= alpha,
                                  "H0 ditolak, tidak homogen",
                                  "Tidak konsisten, periksa"))
  )
}

tabel_homogen <- bind_rows(
  uji_homogen("3.1 Antar tahun",                rfh ~ tahun,      dp),
  uji_homogen("3.2 Jawa vs Luar Jawa",          rfh ~ region,     dp),
  uji_homogen("3.3 Antar 6 pulau",              rfh ~ pulau_6,    dp),
  uji_homogen("3.4 Antar 34 provinsi",          rfh ~ provinsi,   dp),
  uji_homogen("3.5 Kalbar antar tahun (dekadal)", rfh ~ tahun,    kalbar_tahunan),
  uji_homogen("3.5 Kalbar vs Non-Kalbar",       rfh ~ is_kalbar,  dp)
)

kable(tabel_homogen, align = "lccl",
      caption = "Ringkasan Uji Homogenitas Ragam (α = 0.05)")
Ringkasan Uji Homogenitas Ragam (α = 0.05)
Sub-bagian Bartlett p Levene p Keputusan
3.1 Antar tahun 0.105 0.114 H0 diterima, homogen
3.2 Jawa vs Luar Jawa 0.571 0.645 H0 diterima, homogen
3.3 Antar 6 pulau 0.013 0.011 H0 ditolak, tidak homogen
3.4 Antar 34 provinsi 0.900 0.904 H0 diterima, homogen
3.5 Kalbar antar tahun (dekadal) 0.957 0.950 H0 diterima, homogen
3.5 Kalbar vs Non-Kalbar 0.783 0.398 H0 diterima, homogen

4. Uji Stasioneritas Kalbar

adf.test(y)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -2.7837, Lag order = 5, p-value = 0.2492
## alternative hypothesis: stationary
kpss.test(y, null="Level")
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 0.31906, Truncation lag parameter = 4, p-value = 0.1
y_diff <- diff(y)
plot(y_diff, main="Kalbar Differencing", col="darkblue")

adf.test(y_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y_diff
## Dickey-Fuller = -7.7452, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
library(dplyr)
library(tseries)
library(knitr)

adf1  <- adf.test(y)
kpss1 <- kpss.test(y, null = "Level")
adf2  <- adf.test(y_diff)

# format p-value: tseries membatasi p di 0.01 dan 0.10 (nilai batas tabel)
fmt_p <- function(p, lo = 0.01, hi = 0.10) {
  if (p <= lo) "≤ 0.01" else if (p >= hi && hi == 0.10 && p == 0.1) "≥ 0.10" else sprintf("%.3f", p)
}

tabel_stasioner <- tibble(
  Uji = c("ADF (H0: tidak stasioner)",
          "KPSS (H0: stasioner)",
          "ADF setelah diff"),
  Statistik = round(c(adf1$statistic, kpss1$statistic, adf2$statistic), 3),
  `p-value` = c(sprintf("%.3f", adf1$p.value),
                "≥ 0.10",
                "≤ 0.01"),
  Arti = c(
    ifelse(adf1$p.value  < 0.05, "Tolak H0, stasioner", "Gagal tolak H0, mengarah tidak stasioner"),
    ifelse(kpss1$p.value < 0.05, "Tolak H0, tidak stasioner", "Gagal tolak H0, mengarah stasioner"),
    ifelse(adf2$p.value  < 0.05, "Tolak H0, stasioner", "Gagal tolak H0, mengarah tidak stasioner")
  )
)

kable(tabel_stasioner, align = "lccl",
      caption = "Ringkasan Uji Stasioneritas Kalbar (α = 0.05)")
Ringkasan Uji Stasioneritas Kalbar (α = 0.05)
Uji Statistik p-value Arti
ADF (H0: tidak stasioner) -2.784 0.249 Gagal tolak H0, mengarah tidak stasioner
KPSS (H0: stasioner) 0.319 ≥ 0.10 Gagal tolak H0, mengarah stasioner
ADF setelah diff -7.745 ≤ 0.01 Tolak H0, stasioner

5. Analisis Lanjutan

5.1 ANOVA Antar Tahun

summary(aov(rfh ~ tahun, data=dp))
##              Df Sum Sq Mean Sq F value Pr(>F)    
## tahun         4  23113    5778   26.98 <2e-16 ***
## Residuals   165  35339     214                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hasil ANOVA menunjukkan rata-rata curah hujan berbeda nyata antar tahun (F = 26.98; p < 0.001). ## 5.2 ANOVA 6 Pulau + Tukey

model6 <- aov(rfh ~ pulau_6, data=dp)
summary(model6)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## pulau_6       5  16153    3231   12.53 2.64e-10 ***
## Residuals   164  42299     258                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(model6)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rfh ~ pulau_6, data = dp)
## 
## $pulau_6
##                               diff         lwr        upr     p adj
## Jawa-Bali-Nusa           17.366398   2.7211947 32.0116009 0.0101153
## Kalimantan-Bali-Nusa     34.131385  19.0058840 49.2568854 0.0000000
## Maluku-Papua-Bali-Nusa   32.563869  16.7452421 48.3824967 0.0000003
## Sulawesi-Bali-Nusa       12.961903  -1.6833001 27.6071062 0.1153662
## Sumatera-Bali-Nusa       23.087150   9.4532076 36.7210919 0.0000360
## Kalimantan-Jawa          16.764987   4.2235842 29.3063896 0.0022611
## Maluku-Papua-Jawa        15.197472   1.8282914 28.5666518 0.0158573
## Sulawesi-Jawa            -4.404495 -16.3622530  7.5532635 0.8955606
## Sumatera-Jawa             5.720752  -4.9745922 16.4160961 0.6373988
## Maluku-Papua-Kalimantan  -1.567515 -15.4611749 12.3261443 0.9995084
## Sulawesi-Kalimantan     -21.169482 -33.7108843 -8.6280790 0.0000385
## Sumatera-Kalimantan     -11.044235 -22.3883605  0.2998906 0.0612646
## Sulawesi-Maluku-Papua   -19.601966 -32.9711465 -6.2327862 0.0005511
## Sumatera-Maluku-Papua    -9.476720 -21.7297757  2.7763364 0.2296494
## Sumatera-Sulawesi        10.125247  -0.5700974 20.8205908 0.0748245

F = 12.53, p < 0.001, jadi rata-rata berbeda nyata antar pulau (Welch di 3.3 juga sama, p < 0.001). Tukey: Bali-Nusa paling kering, berbeda nyata dari Jawa, Kalimantan, Maluku-Papua, dan Sumatera (kecuali Sulawesi, p = 0.115). Kalimantan dan Maluku-Papua paling basah, keduanya berbeda nyata dari Sulawesi dan Jawa, tetapi tidak berbeda satu sama lain (p = 0.999). Jawa, Sulawesi, dan Sumatera tidak berbeda nyata satu sama lain.

5.3 ANOVA 34 Provinsi (Kruskal-Wallis + Top 5)

model34 <- aov(rfh ~ provinsi, data=dp)
summary(model34)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## provinsi     33  28250   856.1   3.855 1.48e-08 ***
## Residuals   136  30201   222.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dp_34
kruskal.test(rfh ~ provinsi, data=dp) # non-parametrik  untuk 34 provinsi
## 
##  Kruskal-Wallis rank sum test
## 
## data:  rfh by provinsi
## Kruskal-Wallis chi-squared = 79.488, df = 33, p-value = 1.041e-05

ANOVA dan Kruskal-Wallis sama-sama menunjukkan perbedaan curah hujan yang nyata antar 34 provinsi (p < 0.001).

5.4 T-Test Jawa vs Luar Jawa

t.test(rfh ~ region, data=dp, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  rfh by region
## t = -1.1782, df = 168, p-value = 0.2404
## alternative hypothesis: true difference in means between group Jawa and group Luar Jawa is not equal to 0
## 95 percent confidence interval:
##  -11.781532   2.974663
## sample estimates:
##      mean in group Jawa mean in group Luar Jawa 
##                75.28889                79.69232

Tidak terdapat perbedaan rata-rata curah hujan yang nyata antara Jawa dan Luar Jawa (t = -1.178; p = 0.240).

5.5 Forecasting Kalbar

# ARIMA
fit_arima <- auto.arima(y)
summary(fit_arima)
## Series: y 
## ARIMA(1,0,1)(1,0,0)[36] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     sar1      mean
##       0.8660  -0.6180  -0.1233  101.2118
## s.e.  0.0751   0.1058   0.0895    7.3632
## 
## sigma^2 = 1436:  log likelihood = -847.42
## AIC=1704.84   AICc=1705.21   BIC=1720.46
## 
## Training set error measures:
##                      ME    RMSE      MAE       MPE     MAPE      MASE
## Training set -0.1627714 37.4349 29.89787 -27.31146 47.69661 0.6731422
##                    ACF1
## Training set 0.01991346
forecast_arima <- forecast(fit_arima, h=36) 
print(forecast_arima)
##          Point Forecast     Lo 80     Hi 80       Lo 95    Hi 95
## 2026.667       46.12788 -2.428401  94.68417 -28.1325450 120.3883
## 2026.694       52.57044  2.543103 102.59777 -23.9397680 129.0806
## 2026.722       70.34165 19.238944 121.44435  -7.8131939 148.4965
## 2026.750       72.74259 20.848076 124.63710  -6.6232190 152.1084
## 2026.778       66.82680 14.346346 119.30725 -13.4351270 147.0887
## 2026.806       77.45648 24.540883 130.37208  -3.4709416 158.3839
## 2026.833       80.36357 27.123995 133.60315  -1.0593360 161.7865
## 2026.861       81.06697 27.585720 134.54821  -0.7255419 162.8595
## 2026.889       81.50294 27.841174 135.16470  -0.5656465 163.5715
## 2026.917       84.34381 30.547076 138.14054   2.0688058 166.6188
## 2026.944       84.53133 30.633603 138.42905   2.1018702 166.9608
## 2026.972       80.69175 26.718416 134.66509  -1.8533435 163.2368
## 2027.000       87.64370 33.613731 141.67367   5.0119928 170.2754
## 2027.028      101.45299 47.380592 155.52539  18.7563929 184.1496
## 2027.056       99.25126 45.147068 153.35546  16.5060366 181.9965
## 2027.083       96.92839 42.800358 151.05641  14.1467112 179.7101
## 2027.111       89.68027 35.534374 143.82616   6.8712699 172.4893
## 2027.139       98.43906 44.279778 152.59835  15.6095842 181.2685
## 2027.167       95.05518 40.885855 149.22451  12.2103456 177.9000
## 2027.194      100.85846 46.681599 155.03531  18.0021036 183.7148
## 2027.222      100.75330 46.570800 154.93580  17.8883168 183.6183
## 2027.250       96.07551 41.888772 150.26224  13.2040476 178.9470
## 2027.278      101.96631 47.776404 156.15622  19.0899997 184.8426
## 2027.306       99.70330 45.511011 153.89559  16.8233462 182.5833
## 2027.333      100.48954 46.295462 154.68361  17.6068527 183.3722
## 2027.361       97.44444 43.249028 151.63985  14.5597100 180.3292
## 2027.389      100.55425 46.357831 154.75066  17.6679814 183.4405
## 2027.417      105.01044 50.813270 159.20761  22.1230224 187.8979
## 2027.444      102.20824 48.010510 156.40598  19.3199635 185.0965
## 2027.472      105.43281 51.234651 159.63096  22.5438805 188.3217
## 2027.500      110.60325 56.404779 164.80173  27.7138407 193.4927
## 2027.528      110.35672 56.158013 164.55544  27.4669481 193.2465
## 2027.556      108.79001 54.591121 162.98890  25.8999622 191.6801
## 2027.583      111.69732 57.498295 165.89634  28.8070650 194.5876
## 2027.611      111.22386 57.024731 165.42298  28.3334480 194.1143
## 2027.639      110.80785 56.608649 165.00705  27.9173258 193.6984
plot(forecast_arima, main="Forecast Curah Hujan Dekadal Kalbar dengan ARIMA 1 Tahun",
     ylab="Curah Hujan (mm)", xlab="Dekad (36 = 1 tahun)", col="darkgreen")

#  (MAPE kecil = bagus)
accuracy(forecast_arima)
##                      ME    RMSE      MAE       MPE     MAPE      MASE
## Training set -0.1627714 37.4349 29.89787 -27.31146 47.69661 0.6731422
##                    ACF1
## Training set 0.01991346

Model ARIMA(1,0,1)(1,0,0)[36] menghasilkan ramalan yang naik dari sekitar 46 mm menuju rata-rata jangka panjang sekitar 100-110 mm. Interval prediksi yang lebar dan MAPE 47.7% menunjukkan ketidakpastian tinggi, sehingga hasil hanya menggambarkan kecenderungan umum

6. Kesimpulan

  1. Uji homogenitas. Ragam curah hujan homogen antar tahun (Bartlett p = 0.105; Levene p = 0.114), antara Jawa dan Luar Jawa (p = 0.571; 0.645), antar 34 provinsi (p = 0.900; 0.904), serta pada Kalbar antar tahun (p = 0.957; 0.950) dan Kalbar vs Non-Kalbar (p = 0.783; 0.398). Ragam hanya tidak homogen antar 6 pulau (Bartlett p = 0.013; Levene p = 0.011). Kemungkinan penyebabnya, satu pulau memuat provinsi dengan rata-rata yang berbeda-beda.

  2. Stasioneritas Kalbar. Uji ADF (p = 0.249) dan KPSS (p ≥ 0.10) memberi hasil yang bertentangan. Dengan mempertimbangkan bahwa ADF lemah pada data persisten dan auto.arima memilih d = 0, data curah hujan dekadal Kalbar dianggap cenderung stasioner terhadap rata-rata. Differencing tidak digunakan pada model.

  3. Perbandingan rata-rata.

    Rata-rata berbeda nyata antar tahun (F = 26.98; p < 0.001). Nilai 2026 perlu ditafsirkan hati-hati karena data baru sampai Agustus. Rata-rata berbeda nyata antar pulau (Welch p < 0.001). Uji Tukey menunjukkan Bali-Nusa paling kering, Kalimantan dan Maluku-Papua paling basah, sedangkan Jawa, Sulawesi, dan Sumatera tidak berbeda nyata satu sama lain. Rata-rata berbeda nyata antar 34 provinsi (ANOVA p < 0.001; Kruskal-Wallis p < 0.001). Tidak ada perbedaan nyata antara Jawa (75.29 mm) dan Luar Jawa (79.69 mm) (t = -1.178; p = 0.240).

  4. Forecasting Kalbar. Model ARIMA(1,0,1)(1,0,0)[36] meramalkan curah hujan naik dari sekitar 46 mm menuju rata-rata jangka panjang sekitar 100-110 mm. Komponen musimannya tidak signifikan, galatnya cukup besar (RMSE 37.43; MAPE 47.7%), dan interval prediksinya sangat lebar sampai memuat nilai negatif. Hasil ini hanya menggambarkan kecenderungan umum, bukan prediksi akurat per dekad.

7. Referensi