Pendahuluan

Deskripsi Dataset

Dataset yang digunakan adalah data jumlah kejadian bencana alam menurut jenis bencana dan kecamatan di Kota Malang selama periode 2016–2025. Data bersumber dari Badan Pusat Statistik (BPS) Kota Malang dan dapat diakses melalui tautan berikut:

https://malangkota.bps.go.id/id/statistics-table/2/NDY0IzI=/jumlah-bencana-menurut-jenis-bencana-dan-kecamatan-di-kota-malang.html

Dataset terdiri atas 50 observasi yang merupakan kombinasi 5 kecamatan dan 10 tahun pengamatan. Variabel yang tersedia meliputi:

Tahun : tahun pengamatan
Kecamatan : wilayah kecamatan di Kota Malang
Longsor : jumlah kejadian tanah longsor
Cuaca : jumlah kejadian cuaca ekstrem
Banjir : jumlah kejadian banjir
Gempa : jumlah kejadian gempa bumi
Total : total seluruh kejadian bencana

Pada analisis ini, variabel Banjir dipilih sebagai fokus utama karena memiliki variasi kejadian yang cukup besar antar kecamatan dan antar tahun serta merupakan salah satu bencana hidrometeorologis yang berkaitan erat dengan kondisi lingkungan. Analisis dilakukan untuk melihat apakah variasi kejadian banjir antar kecamatan bersifat homogen atau heterogen, menguji apakah terdapat perbedaan rata-rata kejadian banjir berdasarkan kecamatan dan tahun, serta melakukan peramalan jumlah kejadian banjir untuk beberapa tahun berikutnya berdasarkan pola historis tahun 2016–2025.

Tujuan Analisis

Analisis ini bertujuan untuk:

Mendeskripsikan karakteristik jumlah kejadian banjir pada lima kecamatan di Kota Malang selama periode 2016–2025.
Mengidentifikasi homogenitas atau heterogenitas variasi kejadian banjir antar kecamatan melalui visualisasi boxplot dan uji homogenitas varians menggunakan Levene Test berbasis median.
Menguji apakah terdapat perbedaan rata-rata jumlah kejadian banjir antar kecamatan dengan mempertimbangkan variasi antar tahun menggunakan Analysis of Variance (ANOVA).
Mengevaluasi asumsi model ANOVA melalui pemeriksaan normalitas residual dan pola residual.
Mengidentifikasi pasangan kecamatan yang memiliki perbedaan rata-rata kejadian banjir melalui uji lanjut Tukey HSD apabila hasil ANOVA menunjukkan perbedaan yang signifikan.
Menganalisis pola perkembangan total kejadian banjir Kota Malang dari tahun 2016 hingga 2025 sebagai dasar analisis deret waktu.
Menganalisis kecenderungan total kejadian banjir menggunakan Single Moving Average (SMA) orde 3 dan orde 5.
Melakukan peramalan total kejadian banjir Kota Malang periode 2026–2028 menggunakan Double Moving Average (DMA) orde 3.

Persiapan Package

library(readxl)

## Warning: package 'readxl' was built under R version 4.5.3

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(knitr)

Input Data

df <- read_excel("Data_Soal_Statistika_Lingkungan.xlsx")

df <- df %>%
  mutate(
    Tahun = as.integer(Tahun),
    Kecamatan = as.factor(Kecamatan)
  )

# Struktur data
str(df)

## tibble [50 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Tahun        : int [1:50] 2016 2016 2016 2016 2016 2017 2017 2017 2017 2017 ...
##  $ Kecamatan    : Factor w/ 5 levels "Blimbing","Kedungkandang",..: 2 5 3 1 4 2 5 3 1 4 ...
##  $ Tanah Longsor: num [1:50] 8 10 9 5 4 16 24 13 14 10 ...
##  $ Cuaca Ekstrim: num [1:50] 4 1 2 4 1 12 4 15 9 3 ...
##  $ Banjir       : num [1:50] 0 0 1 1 0 0 2 3 4 4 ...
##  $ Gempa Bumi   : num [1:50] 1 1 1 1 1 2 2 2 2 2 ...
##  $ Total Bencana: num [1:50] 13 12 13 11 6 30 32 33 29 19 ...

# Data awal
head(df)

## # A tibble: 6 × 7
##   Tahun Kecamatan     `Tanah Longsor` `Cuaca Ekstrim` Banjir `Gempa Bumi`
##   <int> <fct>                   <dbl>           <dbl>  <dbl>        <dbl>
## 1  2016 Kedungkandang               8               4      0            1
## 2  2016 Sukun                      10               1      0            1
## 3  2016 Klojen                      9               2      1            1
## 4  2016 Blimbing                    5               4      1            1
## 5  2016 Lowokwaru                   4               1      0            1
## 6  2017 Kedungkandang              16              12      0            2
## # ℹ 1 more variable: `Total Bencana` <dbl>

# Dimensi data
dim(df)

## [1] 50  7

# Missing value
colSums(is.na(df))

##         Tahun     Kecamatan Tanah Longsor Cuaca Ekstrim        Banjir 
##             0             0             0             0             0 
##    Gempa Bumi Total Bencana 
##             0             0

Statistik Deskriptif

deskriptif_banjir <- df %>%
  group_by(Kecamatan) %>%
  summarise(
    n = n(),
    Mean = mean(Banjir, na.rm = TRUE),
    Median = median(Banjir, na.rm = TRUE),
    SD = sd(Banjir, na.rm = TRUE),
    Varians = var(Banjir, na.rm = TRUE),
    Min = min(Banjir, na.rm = TRUE),
    Max = max(Banjir, na.rm = TRUE),
    IQR = IQR(Banjir, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    across(
      where(is.numeric),
      ~ round(.x, 2)
    )
  )

knitr::kable(
  deskriptif_banjir,
  caption = "Statistik Deskriptif Kejadian Banjir per Kecamatan"
)

Statistik Deskriptif Kejadian Banjir per Kecamatan
Kecamatan	n	Mean	Median	SD	Varians	Min	Max	IQR
Blimbing	10	15.6	13.5	13.34	178.04	1	35	24.25
Kedungkandang	10	10.3	7.0	11.86	140.68	0	35	17.00
Klojen	10	17.8	9.0	21.41	458.18	1	66	27.75
Lowokwaru	10	17.5	15.5	15.43	238.06	0	52	17.25
Sukun	10	10.6	8.0	12.47	155.60	0	40	14.25

Berdasarkan statistik deskriptif, terdapat perbedaan karakteristik jumlah kejadian banjir pada lima kecamatan di Kota Malang.

Kecamatan Klojen memiliki rata-rata kejadian banjir tertinggi sebesar 17,8 kejadian, diikuti oleh Lowokwaru sebesar 17,5 kejadian dan Blimbing sebesar 15,6 kejadian. Sementara itu, Sukun dan Kedungkandang memiliki rata-rata yang lebih rendah, yaitu masing-masing sebesar 10,6 dan 10,3 kejadian.

Jika dilihat berdasarkan median, Lowokwaru memiliki median tertinggi sebesar 15,5, sedangkan median Klojen sebesar 9. Perbedaan yang cukup besar antara mean dan median pada Klojen mengindikasikan bahwa nilai rata-ratanya dipengaruhi oleh beberapa pengamatan dengan jumlah kejadian yang tinggi. Hal tersebut juga terlihat dari nilai maksimum Klojen yang mencapai 66 kejadian.

Dari sisi penyebaran, Klojen memiliki simpangan baku tertinggi sebesar 21,41, varians sebesar 458,18, dan IQR sebesar 27,75. Hal ini menunjukkan bahwa secara deskriptif Klojen memiliki variasi kejadian banjir yang paling besar dibandingkan kecamatan lainnya.

Kedungkandang memiliki simpangan baku terendah sebesar 11,86. Meskipun secara deskriptif terlihat adanya perbedaan tingkat variasi antar kecamatan, informasi tersebut belum cukup untuk menyimpulkan bahwa varians populasi berbeda secara signifikan. Oleh karena itu, analisis dilanjutkan dengan pengujian homogenitas varians.

Visualisasi

ggplot(
  df,
  aes(
    x = Kecamatan,
    y = Banjir,
    fill = Kecamatan
  )
) +
  geom_boxplot(
    alpha = 0.7,
    outlier.colour = "red"
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.6
  ) +
  labs(
    title = "Sebaran Kejadian Banjir per Kecamatan",
    subtitle = "Kota Malang Tahun 2016-2025",
    x = "Kecamatan",
    y = "Jumlah Kejadian Banjir"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

Boxplot menunjukkan bahwa posisi median dan tingkat penyebaran jumlah kejadian banjir berbeda antar kecamatan. Klojen tampak memiliki penyebaran yang relatif lebih lebar serta memiliki pengamatan dengan jumlah kejadian yang jauh lebih tinggi dibandingkan sebagian besar observasi lainnya. Lowokwaru dan Sukun juga menunjukkan beberapa pengamatan dengan nilai yang relatif tinggi.

Secara visual terdapat indikasi perbedaan tingkat variasi kejadian banjir antar kecamatan. Namun, perbedaan lebar box dan keberadaan nilai ekstrem pada boxplot tidak secara otomatis membuktikan bahwa varians antar kecamatan berbeda secara statistik.

Oleh karena itu, kesimpulan mengenai homogenitas atau heterogenitas varians selanjutnya ditentukan menggunakan Levene Test berbasis median.

Uji Homogenitas

Pengujian homogenitas dilakukan untuk mengetahui apakah varians jumlah kejadian banjir pada lima kecamatan dapat dianggap sama atau terdapat perbedaan varians antar kecamatan.

Hipotesis yang digunakan adalah:

\[H_0:\sigma_1^2=\sigma_2^2=\sigma_3^2=\sigma_4^2=\sigma_5^2\]

\[H_1:\text{minimal terdapat satu varians yang berbeda}\]

Dengan taraf signifikansi:

\[\alpha=0{,}05\]

Daerah penolakan:

\[\text{Tolak }H_0\text{ jika }p\text{-value}<0{,}05\]

uji_homogenitas <- car::leveneTest(
  Banjir ~ Kecamatan,
  data = df
)

uji_homogenitas

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  0.7587 0.5576
##       45

p_levene <- uji_homogenitas$`Pr(>F)`[1]

cat(
  "\nP-value =", round(p_levene, 4), "\n"
)

## 
## P-value = 0.5576

if (p_levene < 0.05) {
  
  cat(
    "Keputusan: Tolak H0\n",
    "Kesimpulan: Varians kejadian banjir antar kecamatan tidak homogen.\n"
  )
  
} else {
  
  cat(
    "Keputusan: Gagal menolak H0\n",
    "Kesimpulan: Belum terdapat bukti adanya perbedaan varians ",
    "kejadian banjir antar kecamatan.\n",
    "Varians antar kecamatan dapat dianggap homogen.\n"
  )
}

## Keputusan: Gagal menolak H0
##  Kesimpulan: Belum terdapat bukti adanya perbedaan varians  kejadian banjir antar kecamatan.
##  Varians antar kecamatan dapat dianggap homogen.

Berdasarkan Levene Test berbasis median diperoleh nilai F sebesar 0,7587 dengan p-value sebesar 0,5576.

Karena p-value sebesar 0,5576 lebih besar dari taraf signifikansi 0,05, maka keputusan yang diperoleh adalah gagal menolak \(H_0\).

Dengan demikian, belum terdapat bukti yang cukup untuk menyatakan bahwa varians jumlah kejadian banjir berbeda antar kecamatan. Oleh karena itu, varians kejadian banjir antar kecamatan dapat dianggap homogen.

Hasil ini menunjukkan bahwa meskipun secara visual terdapat perbedaan tingkat penyebaran pada boxplot, perbedaan tersebut tidak signifikan secara statistik pada taraf signifikansi 5%.

Analisis Anova

ANOVA digunakan untuk menguji apakah terdapat perbedaan rata-rata kejadian banjir antar kecamatan dengan tetap memperhitungkan variasi antar tahun. Model yang digunakan adalah model aditif dengan Kecamatan sebagai faktor utama dan Tahun sebagai faktor waktu.

\[Y_{ij} = \mu + \alpha_i + \beta_j + \varepsilon_{ij}\] Dengan,

\(Y_{ij}\) = jumlah kejadian banjir pada kecamatan ke-i dan tahun ke-j

\(\mu\) = Rata-rata

\(\alpha_i\) = pengaruh kecamatan ke-i

\(\beta_j\) = pengaruh tahun ke-j

\(\varepsilon_{ij}\) = galat/error pengamatan

Model tidak memasukkan interaksi Kecamatan × Tahun. Artinya, analisis berfokus pada pengaruh utama kecamatan dan tahun.

model_anova <- aov(
  Banjir ~ Kecamatan + factor(Tahun),
  data = df
)

summary(model_anova)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Kecamatan      4    539   134.6   3.247   0.0226 *  
## factor(Tahun)  9   9042  1004.7  24.231 1.02e-12 ***
## Residuals     36   1493    41.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pengaruh Kecamatan

\[H_0:\mu_1=\mu_2=\mu_3=\mu_4=\mu_5\]

\[H_1:\text{minimal terdapat satu rata-rata kecamatan yang berbeda}\] Untuk faktor Kecamatan diperoleh F = 3,247 dengan p-value = 0,0226. Karena 0,0226 < 0,05, maka H₀ ditolak. Artinya, terdapat perbedaan rata-rata jumlah kejadian banjir yang signifikan antar kecamatan setelah variasi antar tahun diperhitungkan.

Pengaruh Tahun

Untuk faktor Tahun diperoleh F = 24,231 dengan p-value = \(1,02×10^{-12}\). Nilai tersebut jauh lebih kecil dari 0,05, sehingga terdapat perbedaan kejadian banjir yang signifikan antar tahun selama periode 2016–2025. Temuan ini menunjukkan bahwa variasi temporal merupakan komponen yang sangat kuat dalam data.

UjiNormalitas Residual ANOVA

\[H_0:\text{residual mengikuti distribusi normal}\]

\[H_1:\text{residual tidak mengikuti distribusi normal}\]

residual_anova <- residuals(model_anova)

# Q-Q Plot
qqnorm(
  residual_anova,
  main = "Q-Q Plot Residual ANOVA"
)

qqline(
  residual_anova,
  col = "red",
  lwd = 2
)

# Shapiro-Wilk
uji_normalitas <- shapiro.test(
  residual_anova
)

uji_normalitas

## 
##  Shapiro-Wilk normality test
## 
## data:  residual_anova
## W = 0.95367, p-value = 0.04836

if (uji_normalitas$p.value < 0.05) {
  
  cat(
    "\nKeputusan: Tolak H0\n",
    "Kesimpulan: Residual tidak mengikuti distribusi normal.\n"
  )
  
} else {
  
  cat(
    "\nKeputusan: Gagal menolak H0\n",
    "Kesimpulan: Belum terdapat bukti bahwa residual ",
    "menyimpang dari distribusi normal.\n"
  )
}

## 
## Keputusan: Tolak H0
##  Kesimpulan: Residual tidak mengikuti distribusi normal.

Karena p-value sebesar 0,04836 lebih kecil dari 0,05, maka H₀ ditolak. Secara formal residual model ANOVA belum mengikuti distribusi normal pada taraf signifikansi 5%. Q-Q plot menunjukkan sebagian besar residual di bagian tengah mengikuti garis referensi, tetapi terdapat penyimpangan pada bagian ekor distribusi.

Nilai p-value sangat dekat dengan 0,05, sehingga pelanggaran normalitas perlu dicatat sebagai keterbatasan analisis. Hasil ANOVA tetap dapat disajikan untuk menjawab tugas, tetapi interpretasi inferensialnya dilakukan secara hati-hati dan tidak boleh dianggap tanpa keterbatasan.

RESIDUAL VS FITTED

plot(
  fitted(model_anova),
  residual_anova,
  pch = 16,
  xlab = "Fitted Value",
  ylab = "Residual",
  main = "Residual vs Fitted"
)

abline(
  h = 0,
  col = "red",
  lwd = 2
)

Plot residual terhadap fitted value menunjukkan residual umumnya tersebar di sekitar garis nol, tetapi pada fitted value yang lebih tinggi terlihat penyebaran residual yang lebih besar dan terdapat beberapa observasi yang cukup jauh dari garis nol. Hal ini mengindikasikan model belum menangkap seluruh variasi data secara sempurna.

Meskipun demikian, Levene Test sebelumnya menghasilkan p-value 0,5576 sehingga tidak terdapat bukti statistik mengenai perbedaan varians antar kelompok kecamatan. Oleh karena itu, pemeriksaan visual residual sebaiknya dibaca bersama hasil uji formal, bukan digunakan sendirian untuk menyatakan adanya heteroskedastisitas.

UJI LANJUT TUKEY HSD

Karena faktor Kecamatan signifikan pada ANOVA, dilakukan Tukey HSD untuk mengidentifikasi pasangan kecamatan yang berbeda secara signifikan setelah koreksi pengujian berganda.

p_kecamatan <- summary(model_anova)[[1]][
  "Kecamatan",
  "Pr(>F)"
]

if (
  !is.na(p_kecamatan) &&
  p_kecamatan < 0.05 &&
  p_levene >= 0.05
) {
  
  cat(
    "\nEfek Kecamatan signifikan.",
    "\nDilanjutkan dengan Tukey HSD:\n"
  )
  
  tukey_kecamatan <- TukeyHSD(
    model_anova,
    which = "Kecamatan"
  )
  
  print(tukey_kecamatan)
  
} else {
  
  cat(
    "\nTukey HSD tidak dilakukan karena syarat ",
    "uji lanjut belum terpenuhi.\n"
  )
}

## 
## Efek Kecamatan signifikan. 
## Dilanjutkan dengan Tukey HSD:
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Banjir ~ Kecamatan + factor(Tahun), data = df)
## 
## $Kecamatan
##                         diff        lwr       upr     p adj
## Kedungkandang-Blimbing  -5.3 -13.567129  2.967129 0.3670981
## Klojen-Blimbing          2.2  -6.067129 10.467129 0.9392273
## Lowokwaru-Blimbing       1.9  -6.367129 10.167129 0.9636423
## Sukun-Blimbing          -5.0 -13.267129  3.267129 0.4254774
## Klojen-Kedungkandang     7.5  -0.767129 15.767129 0.0907730
## Lowokwaru-Kedungkandang  7.2  -1.067129 15.467129 0.1130635
## Sukun-Kedungkandang      0.3  -7.967129  8.567129 0.9999720
## Lowokwaru-Klojen        -0.3  -8.567129  7.967129 0.9999720
## Sukun-Klojen            -7.2 -15.467129  1.067129 0.1130635
## Sukun-Lowokwaru         -6.9 -15.167129  1.367129 0.1396577

Seluruh pasangan memiliki adjusted p-value di atas 0,05. Dengan demikian, tidak terdapat pasangan kecamatan tertentu yang menunjukkan perbedaan rata-rata signifikan pada taraf 5% setelah koreksi perbandingan berganda.

Hasil ini tidak bertentangan dengan ANOVA. ANOVA menguji perbedaan secara global antar seluruh kelompok, sedangkan Tukey HSD menguji pasangan satu per satu dengan penyesuaian multiple comparison. Kesimpulannya, terdapat indikasi perbedaan rata-rata secara keseluruhan, tetapi bukti belum cukup kuat untuk menunjuk pasangan kecamatan tertentu yang berbeda secara signifikan.

DATA TIME SERIES BANJIR

banjir_tahun <- df %>%
  group_by(Tahun) %>%
  summarise(
    Total_Banjir = sum(
      Banjir,
      na.rm = TRUE
    ),
    .groups = "drop"
  )

knitr::kable(
  banjir_tahun,
  caption = "Total Kejadian Banjir Kota Malang per Tahun"
)

Total Kejadian Banjir Kota Malang per Tahun
Tahun	Total_Banjir
2016	2
2017	13
2018	9
2019	14
2020	34
2021	92
2022	98
2023	224
2024	115
2025	117

# Membentuk time series
banjir_ts <- ts(
  banjir_tahun$Total_Banjir,
  start = 2016,
  frequency = 1
)

print(banjir_ts)

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]   2  13   9  14  34  92  98 224 115 117

Periode 2016–2019 menunjukkan jumlah kejadian banjir yang relatif rendah. Mulai 2020 terjadi peningkatan, kemudian jumlah kejadian meningkat tajam pada 2021 dan 2022. Puncak kejadian terjadi pada 2023 sebanyak 224 kejadian.

Setelah puncak 2023, total kejadian turun menjadi 115 pada 2024 dan sedikit meningkat menjadi 117 pada 2025. Pola ini tidak menunjukkan kenaikan yang benar-benar konstan. Data lebih tepat menggambarkan kecenderungan meningkat dalam jangka panjang yang disertai fluktuasi besar, terutama lonjakan 2023.

SINGLE MOVING AVERAGE

SMA orde 3

sma3 <- rollmean(
  banjir_ts,
  k = 3,
  align = "right",
  fill = NA
)

SMA orde 5

sma5 <- rollmean(
  banjir_ts,
  k = 5,
  align = "right",
  fill = NA
)

sma3

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]        NA        NA   8.00000  12.00000  19.00000  46.66667  74.66667
##  [8] 138.00000 145.66667 152.00000

sma5

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]    NA    NA    NA    NA  14.4  32.4  49.4  92.4 112.6 129.2

VISUALISASI SINGLE MOVING AVERAGE

ts.plot(
  banjir_ts,
  sma3,
  sma5,
  col = c(
    "black",
    "steelblue",
    "tomato"
  ),
  lwd = 2,
  main = "Single Moving Average Kejadian Banjir",
  xlab = "Tahun",
  ylab = "Total Kejadian Banjir"
)

legend(
  "topleft",
  legend = c(
    "Data Aktual",
    "SMA Orde 3",
    "SMA Orde 5"
  ),
  col = c(
    "black",
    "steelblue",
    "tomato"
  ),
  lty = 1,
  lwd = 2,
  bty = "n"
)

Single Moving Average Orde 3

Single Moving Average orde 3 menggunakan rata-rata tiga observasi terakhir untuk memperhalus fluktuasi data. Nilai yang diperoleh adalah 8; 12; 19; 46,67; 74,67; 138; 145,67; dan 152. Dua periode pertama bernilai NA karena moving average orde 3 belum dapat dihitung sebelum tersedia tiga observasi.

Pada 2023, data aktual mencapai 224 kejadian, sedangkan SMA orde 3 sebesar 138. Perbedaan ini menunjukkan fungsi pemulusan moving average. Setelah lonjakan 2023, SMA masih meningkat menjadi 145,67 pada 2024 dan 152 pada 2025 walaupun data aktual turun ke sekitar 115–117. Hal ini menunjukkan bahwa SMA memiliki keterlambatan dalam merespons perubahan yang sangat tajam.

Single Moving Average Orde 5

SMA orde 5 menghasilkan nilai 14,4; 32,4; 49,4; 92,4; 112,6; dan 129,2. Karena menggunakan lima observasi sekaligus, SMA orde 5 memberikan garis yang lebih halus dibandingkan SMA orde 3.

Pada 2025, SMA orde 3 sebesar 152 sedangkan SMA orde 5 sebesar 129,2. Hal ini menunjukkan bahwa SMA orde 3 lebih responsif terhadap perubahan data terbaru, sedangkan SMA orde 5 melakukan pemulusan yang lebih kuat. Dalam analisis ini, SMA digunakan terutama untuk membaca kecenderungan umum data, bukan sebagai model forecasting akhir.

DOUBLE MOVING AVERAGE

ORDE 3

m <- 3

# Moving Average pertama
Mp <- rollmean(
  banjir_ts,
  k = m,
  align = "right",
  fill = NA
)

# Moving Average kedua
Mpp <- rollmean(
  Mp,
  k = m,
  align = "right",
  fill = NA
)

# Komponen level
a_t <- 2 * Mp - Mpp

# Komponen trend
b_t <- (2 / (m - 1)) * (Mp - Mpp)

# Menampilkan hasil
Mp

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]        NA        NA   8.00000  12.00000  19.00000  46.66667  74.66667
##  [8] 138.00000 145.66667 152.00000

Mpp

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]        NA        NA        NA        NA  13.00000  25.88889  46.77778
##  [8]  86.44444 119.44444 145.22222

a_t

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]        NA        NA        NA        NA  25.00000  67.44444 102.55556
##  [8] 189.55556 171.88889 158.77778

b_t

## Time Series:
## Start = 2016 
## End = 2025 
## Frequency = 1 
##  [1]        NA        NA        NA        NA  6.000000 20.777778 27.888889
##  [8] 51.555556 26.222222  6.777778

Double Moving Average orde 3 dibentuk dari moving average pertama (M′t) dan moving average kedua (M″t). Pada periode terakhir diperoleh M′t = 152 dan M″t = 145,22.

\[a_t = 2M′_t − M''_t = 157,8\]

\[b_t = [2/(m−1)](M′_t − M''_t) = 6,78\] Nilai \(b_t\) yang positif menunjukkan bahwa berdasarkan pola hasil pemulusan terakhir, model menangkap kecenderungan peningkatan sekitar 6,78 kejadian untuk setiap periode forecast berikutnya.

\[\hat{Y}_{t+k} = a_t + b_t k\] \[\hat{Y}_{t+k} = 158,78 + 6,78k\]

NILAI LEVEL DAN TREND TERAKHIR

a_last <- as.numeric(
  tail(
    na.omit(a_t),
    1
  )
)

b_last <- as.numeric(
  tail(
    na.omit(b_t),
    1
  )
)

cat(
  "Nilai level terakhir (a) =",
  round(a_last, 2),
  "\n"
)

## Nilai level terakhir (a) = 158.78

cat(
  "Nilai trend terakhir (b) =",
  round(b_last, 2),
  "\n"
)

## Nilai trend terakhir (b) = 6.78

FORECAST 2026-2028

DOUBLE MOVING AVERAGE ORDE 3

forecast_2026 <- a_last + b_last * 1
forecast_2027 <- a_last + b_last * 2
forecast_2028 <- a_last + b_last * 3

hasil_forecast <- data.frame(
  Tahun = 2026:2028,
  Forecast = c(
    forecast_2026,
    forecast_2027,
    forecast_2028
  )
)

hasil_forecast <- hasil_forecast %>%
  mutate(
    Forecast = round(
      Forecast,
      2
    )
  )

knitr::kable(
  hasil_forecast,
  caption = "Hasil Forecast Total Kejadian Banjir Tahun 2026-2028 Menggunakan Double Moving Average Orde 3"
)

Hasil Forecast Total Kejadian Banjir Tahun 2026-2028 Menggunakan Double Moving Average Orde 3
Tahun	Forecast
2026	165.56
2027	172.33
2028	179.11

PERBANDINGAN DATA AKTUAL DAN FORECAST

tahun_aktual <- banjir_tahun$Tahun
nilai_aktual <- banjir_tahun$Total_Banjir

tahun_forecast <- hasil_forecast$Tahun
nilai_forecast <- hasil_forecast$Forecast


plot(
  tahun_aktual,
  nilai_aktual,
  type = "o",
  pch = 16,
  lwd = 2,
  xlim = c(2016, 2028),
  ylim = range(
    c(
      nilai_aktual,
      nilai_forecast
    ),
    na.rm = TRUE
  ),
  xlab = "Tahun",
  ylab = "Total Kejadian Banjir",
  main = "Perbandingan Data Aktual dan Forecast Kejadian Banjir"
)

# Garis forecast dimulai dari data aktual terakhir
lines(
  x = c(
    2025,
    tahun_forecast
  ),
  y = c(
    tail(nilai_aktual, 1),
    nilai_forecast
  ),
  type = "o",
  pch = 17,
  lty = 2,
  lwd = 2
)

legend(
  "topleft",
  legend = c(
    "Data Aktual",
    "Forecast DMA Orde 3"
  ),
  lty = c(1, 2),
  pch = c(16, 17),
  lwd = 2,
  bty = "n"
)

Berdasarkan Double Moving Average orde 3, total kejadian banjir diproyeksikan sebesar 165,56 pada 2026, 172,33 pada 2027, dan 179,11 pada 2028. Secara praktis, angka tersebut dapat dibaca sebagai sekitar 166, 172, dan 179 kejadian.

Forecast menunjukkan pola peningkatan bertahap. Dari 2026 hingga 2028, model memperkirakan kenaikan sekitar 13,55 kejadian. Namun, angka forecast bukan merupakan jumlah kejadian yang pasti akan terjadi. Model sepenuhnya dibentuk dari pola historis 2016–2025 dan sangat dipengaruhi oleh struktur data beberapa tahun terakhir, termasuk lonjakan besar pada 2023.

Karena data hanya terdiri atas 10 observasi tahunan, hasil forecasting sebaiknya dipandang sebagai proyeksi berdasarkan kecenderungan historis, bukan prediksi yang bersifat pasti. Interpretasi ini penting agar hasil tidak dilebih-lebihkan melampaui informasi yang didukung oleh data.

PLOT DATA AKTUAL + SMA + FORECAST

plot(
  tahun_aktual,
  nilai_aktual,
  type = "o",
  pch = 16,
  lwd = 2,
  xlim = c(2016, 2028),
  ylim = range(
    c(
      nilai_aktual,
      as.numeric(sma3),
      as.numeric(sma5),
      nilai_forecast
    ),
    na.rm = TRUE
  ),
  xlab = "Tahun",
  ylab = "Total Kejadian Banjir",
  main = "Data Aktual, Moving Average, dan Forecast Banjir"
)

# SMA orde 3
lines(
  x = as.numeric(time(sma3)),
  y = as.numeric(sma3),
  lwd = 2,
  lty = 2
)

# SMA orde 5
lines(
  x = as.numeric(time(sma5)),
  y = as.numeric(sma5),
  lwd = 2,
  lty = 3
)

# Forecast DMA orde 3
lines(
  x = c(
    2025,
    tahun_forecast
  ),
  y = c(
    tail(nilai_aktual, 1),
    nilai_forecast
  ),
  type = "o",
  pch = 17,
  lwd = 2,
  lty = 4
)

legend(
  "topleft",
  legend = c(
    "Data Aktual",
    "SMA Orde 3",
    "SMA Orde 5",
    "Forecast DMA Orde 3"
  ),
  lty = c(
    1,
    2,
    3,
    4
  ),
  pch = c(
    16,
    NA,
    NA,
    17
  ),
  lwd = 2,
  bty = "n"
)

Data aktual 2016–2025 menunjukkan fluktuasi yang cukup besar, terutama lonjakan menjadi 224 kejadian pada 2023 dan penurunan kembali ke 115–117 pada 2024–2025. Sebaliknya, garis forecast DMA orde 3 untuk 2026–2028 lebih halus dan menunjukkan kenaikan bertahap.

Perbedaan karakter tersebut terjadi karena DMA memanfaatkan hasil pemulusan dan tren, sehingga tidak mengikuti setiap fluktuasi ekstrem data aktual. Dengan demikian, grafik gabungan data aktual, SMA, dan forecast membantu menunjukkan hubungan antara data historis, proses pemulusan, dan proyeksi periode mendatang secara visual.

Analisis Homogenitas, ANOVA, dan Forecasting Kejadian Banjir Kota Malang Tahun 2016–2025

Ridla Iasa Ramadhania

2026-09-25