Pemodelan Data Lingkungan Menggunakan GLM dan Regresi Nonlinear
Rahayu Puspita Sari (3338240049)
KASUS 1. PEMODELAN JUMLAH LAPORAN KEJADIAN BANJIR MENGGUNAKAN GENERALIZED LINEAR MODEL
Data dan Sumber Data
Data yang digunakan dalam analisis ini merupakan data laporan kejadian banjir periode Januari 2020 sampai Desember 2025 yang bersumber dari https://data.mendeley.com/datasets/d7j3s64vvr/4.
Data mencakup 27 kabupaten/kota di Provinsi Jawa Barat, yang terdiri dari 18 kabupaten dan 9 kota.
Setiap baris data merepresentasikan satu laporan kejadian banjir yang dilengkapi dengan informasi wilayah dan tanggal kejadian. Dalam analisis ini, data dikelompokkan berdasarkan kabupaten/kota dan bulan untuk memperoleh jumlah laporan kejadian banjir pada setiap wilayah selama periode pengamatan.
Variabel respons yang digunakan adalah jumlah laporan kejadian banjir
(jumlah_banjir), sedangkan variabel prediktor meliputi tren
tahun (tahun_tren), bulan (bulan), dan
kabupaten/kota (district_city).
Persiapan Data
# Membaca data
banjir <- read_excel("C:/Users/rahay/Downloads/Flood Event and Daily Rainfall Data Records across/Flood Event and Daily Rainfall Data Records across/Flood_report_Raw_2020-2025.xlsx")
# Melihat struktur data
str(banjir)## tibble [4,224 × 11] (S3: tbl_df/tbl/data.frame)
## $ id : chr [1:4224] "cf3f0f97-ddb2-490e-b8ee-5e534973591e" "1748f982-d721-4610-ae97-31482fe41fe0" "3d6f2b8c-5628-429b-b994-16f337b0b148" "693454f8-6a66-415d-af21-00cbf3764099" ...
## $ district_city : chr [1:4224] "Sukabumi" "Sukabumi" "Sukabumi" "Sukabumi" ...
## $ sub_district : chr [1:4224] "pabuaran" "Cibadak" "Tegalbuleud" "Purabaya" ...
## $ event_date : chr [1:4224] "2020-02-04T00:00:00" "2020-05-12T00:00:00" "2020-06-22T00:00:00" "2020-06-24T00:00:00" ...
## $ latitude : num [1:4224] -7.16 -6.88 -7.41 -7.1 -7.09 ...
## $ longitude : num [1:4224] 107 107 107 107 107 ...
## $ fatalities : num [1:4224] 0 0 0 0 0 0 0 3 0 0 ...
## $ missing_people : num [1:4224] 0 0 0 0 0 0 0 0 0 0 ...
## $ severely_damaged_house: num [1:4224] 0 0 0 0 0 0 0 84 0 0 ...
## $ damaged_building : num [1:4224] 0 0 0 0 0 0 0 0 0 0 ...
## $ submerged_house : num [1:4224] 0 0 0 0 0 2 1 327 0 0 ...# Melihat beberapa baris pertama
head(banjir)## # A tibble: 6 × 11
## id district_city sub_district event_date latitude longitude fatalities
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 cf3f0f97-… Sukabumi pabuaran 2020-02-0… -7.16 107. 0
## 2 1748f982-… Sukabumi Cibadak 2020-05-1… -6.88 107. 0
## 3 3d6f2b8c-… Sukabumi Tegalbuleud 2020-06-2… -7.41 107. 0
## 4 693454f8-… Sukabumi Purabaya 2020-06-2… -7.10 107. 0
## 5 13e18e30-… Sukabumi jampang ten… 2020-06-2… -7.09 107. 0
## 6 21414c79-… Sukabumi Cisolok 2020-09-0… -6.96 106. 0
## # ℹ 4 more variables: missing_people <dbl>, severely_damaged_house <dbl>,
## # damaged_building <dbl>, submerged_house <dbl># Mengecek missing value
colSums(is.na(banjir))## id district_city sub_district
## 0 0 0
## event_date latitude longitude
## 0 0 0
## fatalities missing_people severely_damaged_house
## 0 0 0
## damaged_building submerged_house
## 0 0# Mengubah format tanggal
banjir$event_date <- as.Date(
substr(as.character(banjir$event_date), 1, 10)
)
# Membuat variabel tahun dan bulan
banjir <- banjir %>%
mutate(
tahun = year(event_date),
bulan = month(event_date)
)
# Menghitung jumlah laporan banjir
data_bulanan <- banjir %>%
count(
district_city,
tahun,
bulan,
name = "jumlah_banjir"
)
# Membuat kalender bulanan
kalender <- expand_grid(
district_city = sort(unique(banjir$district_city)),
tanggal = seq.Date(
from = as.Date("2020-01-01"),
to = as.Date("2025-12-01"),
by = "month"
)
) %>%
mutate(
tahun = year(tanggal),
bulan = month(tanggal),
hari = days_in_month(tanggal)
)
# Menggabungkan data dengan kalender
data_glm <- kalender %>%
left_join(
data_bulanan,
by = c("district_city", "tahun", "bulan")
) %>%
mutate(
jumlah_banjir = replace_na(jumlah_banjir, 0),
tahun_tren = tahun - 2020,
bulan = factor(bulan),
district_city = factor(district_city)
)
# Melihat hasil data
head(data_glm)## # A tibble: 6 × 7
## district_city tanggal tahun bulan hari jumlah_banjir tahun_tren
## <fct> <date> <dbl> <fct> <int> <int> <dbl>
## 1 Bandung 2020-01-01 2020 1 31 3 0
## 2 Bandung 2020-02-01 2020 2 29 6 0
## 3 Bandung 2020-03-01 2020 3 31 3 0
## 4 Bandung 2020-04-01 2020 4 30 1 0
## 5 Bandung 2020-05-01 2020 5 31 1 0
## 6 Bandung 2020-06-01 2020 6 30 0 0# Melihat dimensi data
dim(data_glm)## [1] 1944 7Data awal terdiri dari 4.224 laporan kejadian banjir dengan 11 variabel yang mencakup informasi lokasi, tanggal kejadian, koordinat, korban, dan kerusakan bangunan. Data mencakup 27 kabupaten/kota di Provinsi Jawa Barat selama periode 2020–2025. Hasil pemeriksaan menunjukkan tidak terdapat missing value.
Selanjutnya, data dikelompokkan berdasarkan kabupaten/kota dan bulan untuk memperoleh jumlah laporan banjir. Bulan tanpa laporan diberi nilai 0. Setelah proses pengolahan, diperoleh 1.944 pengamatan dengan 7 variabel yang siap digunakan dalam pemodelan GLM
Statistik Deskriptif
# Statistik deskriptif jumlah banjir
summary(data_glm$jumlah_banjir)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 2.173 2.000 157.000# Rata-rata
mean(data_glm$jumlah_banjir)## [1] 2.17284# Varians
var(data_glm$jumlah_banjir)## [1] 42.29641# Standar deviasi
sd(data_glm$jumlah_banjir)## [1] 6.503569# Jumlah data nol
sum(data_glm$jumlah_banjir == 0)## [1] 1224# Jumlah total laporan banjir
sum(data_glm$jumlah_banjir)## [1] 4224Berdasarkan statistik deskriptif, terdapat 4.224 laporan banjir dengan rata-rata 2,17 laporan per kabupaten/kota setiap bulan. Varians sebesar 42,30 lebih besar dibandingkan rata-ratanya, yang mengindikasikan adanya overdispersion. Selain itu, terdapat 1.224 pengamatan dengan jumlah laporan nol.
# Jumlah laporan berdasarkan tahun
data_glm %>%
group_by(tahun) %>%
summarise(
total_banjir = sum(jumlah_banjir)
)## # A tibble: 6 × 2
## tahun total_banjir
## <dbl> <int>
## 1 2020 427
## 2 2021 477
## 3 2022 710
## 4 2023 700
## 5 2024 794
## 6 2025 1116# Visualisasi Jumlah Laporan Banjir Berdasarkan Tahun
data_glm %>%
group_by(tahun) %>%
summarise(
total_banjir = sum(jumlah_banjir)
) %>%
ggplot(aes(x = tahun, y = total_banjir)) +
geom_col(fill = "steelblue") +
labs(
title = "Jumlah Laporan Kejadian Banjir Berdasarkan Tahun",
x = "Tahun",
y = "Jumlah Laporan Banjir"
) +
theme_minimal()
Berdasarkan grafik, jumlah laporan kejadian banjir cenderung meningkat selama periode 2020–2025. Jumlah laporan terendah terjadi pada tahun 2020 sebanyak 427 laporan, sedangkan tertinggi pada tahun 2025 sebanyak 1.116 laporan. Hal ini menunjukkan adanya peningkatan jumlah laporan banjir dari tahun ke tahun.
# Jumlah laporan berdasarkan bulan
data_glm %>%
group_by(bulan) %>%
summarise(
total_banjir = sum(jumlah_banjir)
)## # A tibble: 12 × 2
## bulan total_banjir
## <fct> <int>
## 1 1 456
## 2 2 672
## 3 3 653
## 4 4 287
## 5 5 298
## 6 6 86
## 7 7 251
## 8 8 91
## 9 9 109
## 10 10 307
## 11 11 457
## 12 12 557# Visualisasi Jumlah Laporan Banjir Berdasarkan Bulan
data_glm %>%
group_by(bulan) %>%
summarise(
total_banjir = sum(jumlah_banjir)
) %>%
ggplot(aes(x = bulan, y = total_banjir)) +
geom_col(fill = "steelblue") +
labs(
title = "Jumlah Laporan Kejadian Banjir Berdasarkan Bulan",
x = "Bulan",
y = "Jumlah Laporan Banjir"
) +
theme_minimal()
Berdasarkan grafik, jumlah laporan banjir tertinggi terjadi pada bulan Februari sebanyak 672 laporan, diikuti Maret sebanyak 653 laporan. Sementara itu, jumlah laporan terendah terjadi pada bulan Juni sebanyak 86 laporan. Secara umum, laporan banjir cenderung lebih banyak pada awal dan akhir tahun, serta lebih sedikit pada pertengahan tahun.
Pemodelan GLM
# MODEL POISSON
model_poisson <- glm(
jumlah_banjir ~ tahun_tren + bulan +
district_city + offset(log(hari)),
family = poisson(link = "log"),
data = data_glm
)
summary(model_poisson)##
## Call:
## glm(formula = jumlah_banjir ~ tahun_tren + bulan + district_city +
## offset(log(hari)), family = poisson(link = "log"), data = data_glm)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.757563 0.069127 -25.425 < 2e-16 ***
## tahun_tren 0.181799 0.009282 19.587 < 2e-16 ***
## bulan2 0.478552 0.060672 7.888 3.08e-15 ***
## bulan3 0.359084 0.061028 5.884 4.01e-09 ***
## bulan4 -0.430221 0.075348 -5.710 1.13e-08 ***
## bulan5 -0.425399 0.074490 -5.711 1.12e-08 ***
## bulan6 -1.635356 0.117562 -13.911 < 2e-16 ***
## bulan7 -0.597040 0.078594 -7.596 3.04e-14 ***
## bulan8 -1.611633 0.114813 -14.037 < 2e-16 ***
## bulan9 -1.398355 0.106617 -13.116 < 2e-16 ***
## bulan10 -0.395645 0.073826 -5.359 8.36e-08 ***
## bulan11 0.034980 0.066190 0.528 0.597165
## bulan12 0.200072 0.063153 3.168 0.001535 **
## district_cityBandung_City -2.118907 0.137789 -15.378 < 2e-16 ***
## district_cityBanjar_City -5.503297 0.708545 -7.767 8.03e-15 ***
## district_cityBekasi -0.304800 0.069275 -4.400 1.08e-05 ***
## district_cityBekasi_City -0.883238 0.083443 -10.585 < 2e-16 ***
## district_cityBogor -0.258908 0.068375 -3.787 0.000153 ***
## district_cityBogor_City -1.707808 0.115207 -14.824 < 2e-16 ***
## district_cityCiamis -1.291169 0.097181 -13.286 < 2e-16 ***
## district_cityCianjur -1.060646 0.088989 -11.919 < 2e-16 ***
## district_cityCimahi_City -1.801995 0.119926 -15.026 < 2e-16 ***
## district_cityCirebon -0.318708 0.069554 -4.582 4.60e-06 ***
## district_cityCirebon_City -1.789725 0.119297 -15.002 < 2e-16 ***
## district_cityDepok_City -1.621733 0.111112 -14.595 < 2e-16 ***
## district_cityGarut -1.020294 0.087672 -11.638 < 2e-16 ***
## district_cityIndramayu -1.814417 0.120568 -15.049 < 2e-16 ***
## district_cityKarawang 0.057385 0.062926 0.912 0.361804
## district_cityKuningan -2.226152 0.144584 -15.397 < 2e-16 ***
## district_cityMajalengka -1.707808 0.115207 -14.824 < 2e-16 ***
## district_cityPangandaran -2.085570 0.135757 -15.362 < 2e-16 ***
## district_cityPurwakarta -2.864240 0.194296 -14.742 < 2e-16 ***
## district_citySubang -1.542484 0.107520 -14.346 < 2e-16 ***
## district_citySukabumi -0.561655 0.074887 -7.500 6.38e-14 ***
## district_citySukabumi_City -1.865711 0.123266 -15.136 < 2e-16 ***
## district_citySumedang -1.653149 0.112583 -14.684 < 2e-16 ***
## district_cityTasikmalaya -2.585526 0.170481 -15.166 < 2e-16 ***
## district_cityTasikmalaya_City -2.085570 0.135757 -15.362 < 2e-16 ***
## district_cityWest_Bandung -2.171092 0.141045 -15.393 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 13932.9 on 1943 degrees of freedom
## Residual deviance: 9054.5 on 1905 degrees of freedom
## AIC: 11348
##
## Number of Fisher Scoring iterations: 7# Rasio deviance
deviance(model_poisson) /
df.residual(model_poisson)## [1] 4.753025# Rasio Pearson
sum(residuals(
model_poisson,
type = "pearson"
)^2) / df.residual(model_poisson)## [1] 9.056242Hasil model menunjukkan bahwa variabel tahun, bulan, dan sebagian besar wilayah berpengaruh signifikan terhadap jumlah laporan banjir. Namun, rasio deviance sebesar 4,75 dan rasio Pearson sebesar 9,06 menunjukkan adanya overdispersion, sehingga model Poisson kurang sesuai untuk data ini.
# MODEL QUASI-POISSON
model_quasi <- glm(
jumlah_banjir ~ tahun_tren + bulan +
district_city + offset(log(hari)),
family = quasipoisson(link = "log"),
data = data_glm
)
summary(model_quasi)##
## Call:
## glm(formula = jumlah_banjir ~ tahun_tren + bulan + district_city +
## offset(log(hari)), family = quasipoisson(link = "log"), data = data_glm)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.75756 0.20803 -8.449 < 2e-16 ***
## tahun_tren 0.18180 0.02793 6.509 9.68e-11 ***
## bulan2 0.47855 0.18258 2.621 0.008837 **
## bulan3 0.35908 0.18365 1.955 0.050703 .
## bulan4 -0.43022 0.22675 -1.897 0.057934 .
## bulan5 -0.42540 0.22417 -1.898 0.057887 .
## bulan6 -1.63536 0.35379 -4.622 4.05e-06 ***
## bulan7 -0.59704 0.23652 -2.524 0.011674 *
## bulan8 -1.61163 0.34551 -4.664 3.31e-06 ***
## bulan9 -1.39836 0.32085 -4.358 1.38e-05 ***
## bulan10 -0.39565 0.22217 -1.781 0.075100 .
## bulan11 0.03498 0.19919 0.176 0.860617
## bulan12 0.20007 0.19005 1.053 0.292597
## district_cityBandung_City -2.11891 0.41466 -5.110 3.54e-07 ***
## district_cityBanjar_City -5.50330 2.13227 -2.581 0.009927 **
## district_cityBekasi -0.30480 0.20847 -1.462 0.143893
## district_cityBekasi_City -0.88324 0.25111 -3.517 0.000446 ***
## district_cityBogor -0.25891 0.20577 -1.258 0.208450
## district_cityBogor_City -1.70781 0.34670 -4.926 9.12e-07 ***
## district_cityCiamis -1.29117 0.29245 -4.415 1.07e-05 ***
## district_cityCianjur -1.06065 0.26780 -3.961 7.75e-05 ***
## district_cityCimahi_City -1.80199 0.36090 -4.993 6.49e-07 ***
## district_cityCirebon -0.31871 0.20931 -1.523 0.128016
## district_cityCirebon_City -1.78972 0.35901 -4.985 6.75e-07 ***
## district_cityDepok_City -1.62173 0.33438 -4.850 1.33e-06 ***
## district_cityGarut -1.02029 0.26384 -3.867 0.000114 ***
## district_cityIndramayu -1.81442 0.36283 -5.001 6.24e-07 ***
## district_cityKarawang 0.05738 0.18937 0.303 0.761898
## district_cityKuningan -2.22615 0.43511 -5.116 3.43e-07 ***
## district_cityMajalengka -1.70781 0.34670 -4.926 9.12e-07 ***
## district_cityPangandaran -2.08557 0.40854 -5.105 3.64e-07 ***
## district_cityPurwakarta -2.86424 0.58471 -4.899 1.05e-06 ***
## district_citySubang -1.54248 0.32357 -4.767 2.01e-06 ***
## district_citySukabumi -0.56165 0.22536 -2.492 0.012779 *
## district_citySukabumi_City -1.86571 0.37095 -5.030 5.38e-07 ***
## district_citySumedang -1.65315 0.33880 -4.879 1.15e-06 ***
## district_cityTasikmalaya -2.58553 0.51304 -5.040 5.11e-07 ***
## district_cityTasikmalaya_City -2.08557 0.40854 -5.105 3.64e-07 ***
## district_cityWest_Bandung -2.17109 0.42446 -5.115 3.45e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 9.056243)
##
## Null deviance: 13932.9 on 1943 degrees of freedom
## Residual deviance: 9054.5 on 1905 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7# Nilai dispersi
summary(model_quasi)$dispersion## [1] 9.056243Model Quasi-Poisson menghasilkan nilai dispersi sebesar 9,06. Model ini mengakomodasi overdispersion dengan menyesuaikan standar error. Hasilnya, beberapa variabel yang sebelumnya signifikan pada model Poisson menjadi tidak signifikan.
# MODEL NEGATIVE BINOMIAL
model_nb <- glm.nb(
jumlah_banjir ~ tahun_tren + bulan +
district_city + offset(log(hari)),
data = data_glm
)
summary(model_nb)##
## Call:
## glm.nb(formula = jumlah_banjir ~ tahun_tren + bulan + district_city +
## offset(log(hari)), data = data_glm, init.theta = 0.3006772968,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.74425 0.27899 -6.252 4.05e-10 ***
## tahun_tren 0.18151 0.02832 6.409 1.46e-10 ***
## bulan2 0.54194 0.22098 2.452 0.01419 *
## bulan3 0.24728 0.22225 1.113 0.26587
## bulan4 -0.39232 0.22834 -1.718 0.08576 .
## bulan5 -0.31703 0.22712 -1.396 0.16275
## bulan6 -1.60538 0.25216 -6.367 1.93e-10 ***
## bulan7 -0.65002 0.23124 -2.811 0.00494 **
## bulan8 -2.01333 0.26503 -7.597 3.04e-14 ***
## bulan9 -1.24472 0.24267 -5.129 2.91e-07 ***
## bulan10 -0.14766 0.22541 -0.655 0.51242
## bulan11 0.10932 0.22349 0.489 0.62475
## bulan12 0.14617 0.22296 0.656 0.51210
## district_cityBandung_City -2.20998 0.34825 -6.346 2.21e-10 ***
## district_cityBanjar_City -5.54277 0.78586 -7.053 1.75e-12 ***
## district_cityBekasi -0.64940 0.31982 -2.031 0.04230 *
## district_cityBekasi_City -1.03372 0.32395 -3.191 0.00142 **
## district_cityBogor 0.16592 0.31445 0.528 0.59775
## district_cityBogor_City -1.50727 0.33117 -4.551 5.33e-06 ***
## district_cityCiamis -1.29266 0.32756 -3.946 7.93e-05 ***
## district_cityCianjur -1.05413 0.32421 -3.251 0.00115 **
## district_cityCimahi_City -1.77508 0.33661 -5.273 1.34e-07 ***
## district_cityCirebon -0.63679 0.31970 -1.992 0.04639 *
## district_cityCirebon_City -1.90845 0.33978 -5.617 1.95e-08 ***
## district_cityDepok_City -1.60032 0.33293 -4.807 1.53e-06 ***
## district_cityGarut -0.71798 0.32046 -2.240 0.02506 *
## district_cityIndramayu -1.95891 0.34106 -5.744 9.27e-09 ***
## district_cityKarawang -0.29503 0.31702 -0.931 0.35203
## district_cityKuningan -2.34155 0.35261 -6.641 3.12e-11 ***
## district_cityMajalengka -1.93548 0.34046 -5.685 1.31e-08 ***
## district_cityPangandaran -1.77705 0.33666 -5.279 1.30e-07 ***
## district_cityPurwakarta -2.80037 0.37172 -7.533 4.94e-14 ***
## district_citySubang -1.83027 0.33788 -5.417 6.06e-08 ***
## district_citySukabumi -0.39542 0.31773 -1.245 0.21331
## district_citySukabumi_City -2.01096 0.34244 -5.872 4.29e-09 ***
## district_citySumedang -1.71627 0.33532 -5.118 3.08e-07 ***
## district_cityTasikmalaya -1.95322 0.34091 -5.729 1.01e-08 ***
## district_cityTasikmalaya_City -1.79990 0.33718 -5.338 9.39e-08 ***
## district_cityWest_Bandung -2.32098 0.35190 -6.596 4.24e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3007) family taken to be 1)
##
## Null deviance: 2077.6 on 1943 degrees of freedom
## Residual deviance: 1470.4 on 1905 degrees of freedom
## AIC: 5854
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3007
## Std. Err.: 0.0166
##
## 2 x log-likelihood: -5773.9560# Rasio Pearson
sum(residuals(
model_nb,
type = "pearson"
)^2) / df.residual(model_nb)## [1] 1.310911Model Negative Binomial menghasilkan nilai AIC sebesar 5.854, lebih rendah dibandingkan Poisson sebesar 11.348. Rasio Pearson juga turun menjadi 1,31, yang menunjukkan bahwa model ini lebih baik dalam menangani overdispersion.
Perbandingan Model
# Perbandingan AIC
AIC(
model_poisson,
model_nb
)## df AIC
## model_poisson 39 11348.300
## model_nb 40 5853.956# Membandingkan nilai dispersi
dispersi <- data.frame(
Model = c(
"Poisson",
"Quasi-Poisson",
"Negative Binomial"
),
Rasio_Pearson = c(
sum(residuals(model_poisson,
type = "pearson")^2) /
df.residual(model_poisson),
summary(model_quasi)$dispersion,
sum(residuals(model_nb,
type = "pearson")^2) /
df.residual(model_nb)
)
)
dispersi## Model Rasio_Pearson
## 1 Poisson 9.056242
## 2 Quasi-Poisson 9.056243
## 3 Negative Binomial 1.310911Berdasarkan nilai AIC, model Negative Binomial memiliki nilai AIC sebesar 5.853,956, lebih rendah dibandingkan model Poisson sebesar 11.348,300. Hal ini menunjukkan bahwa model Negative Binomial lebih baik berdasarkan kriteria AIC.
Berdasarkan rasio Pearson, model Poisson dan Quasi-Poisson memiliki nilai sekitar 9,06, sedangkan Negative Binomial sebesar 1,31. Nilai tersebut menunjukkan bahwa Negative Binomial lebih mampu mengatasi masalah overdispersion pada data.
Kesimpulan: Model Negative Binomial dipilih sebagai model terbaik karena memiliki nilai AIC paling rendah dan rasio Pearson yang lebih mendekati 1.
Interpretasi IRR
# Menghitung IRR
IRR <- exp(coef(model_nb))
# Interval kepercayaan 95%
CI <- exp(confint.default(model_nb))
# Membuat tabel hasil
hasil_irr <- data.frame(
Koefisien = coef(model_nb),
IRR = IRR,
Batas_Bawah = CI[, 1],
Batas_Atas = CI[, 2]
)
round(hasil_irr, 4)## Koefisien IRR Batas_Bawah Batas_Atas
## (Intercept) -1.7442 0.1748 0.1012 0.3020
## tahun_tren 0.1815 1.1990 1.1343 1.2675
## bulan2 0.5419 1.7193 1.1150 2.6513
## bulan3 0.2473 1.2805 0.8283 1.9796
## bulan4 -0.3923 0.6755 0.4318 1.0568
## bulan5 -0.3170 0.7283 0.4667 1.1367
## bulan6 -1.6054 0.2008 0.1225 0.3292
## bulan7 -0.6500 0.5220 0.3318 0.8214
## bulan8 -2.0133 0.1335 0.0794 0.2245
## bulan9 -1.2447 0.2880 0.1790 0.4634
## bulan10 -0.1477 0.8627 0.5546 1.3420
## bulan11 0.1093 1.1155 0.7198 1.7287
## bulan12 0.1462 1.1574 0.7476 1.7917
## district_cityBandung_City -2.2100 0.1097 0.0554 0.2171
## district_cityBanjar_City -5.5428 0.0039 0.0008 0.0183
## district_cityBekasi -0.6494 0.5224 0.2791 0.9777
## district_cityBekasi_City -1.0337 0.3557 0.1885 0.6711
## district_cityBogor 0.1659 1.1805 0.6374 2.1864
## district_cityBogor_City -1.5073 0.2215 0.1157 0.4239
## district_cityCiamis -1.2927 0.2745 0.1445 0.5217
## district_cityCianjur -1.0541 0.3485 0.1846 0.6579
## district_cityCimahi_City -1.7751 0.1695 0.0876 0.3278
## district_cityCirebon -0.6368 0.5290 0.2827 0.9899
## district_cityCirebon_City -1.9085 0.1483 0.0762 0.2887
## district_cityDepok_City -1.6003 0.2018 0.1051 0.3876
## district_cityGarut -0.7180 0.4877 0.2603 0.9140
## district_cityIndramayu -1.9589 0.1410 0.0723 0.2752
## district_cityKarawang -0.2950 0.7445 0.4000 1.3858
## district_cityKuningan -2.3415 0.0962 0.0482 0.1920
## district_cityMajalengka -1.9355 0.1444 0.0741 0.2813
## district_cityPangandaran -1.7771 0.1691 0.0874 0.3272
## district_cityPurwakarta -2.8004 0.0608 0.0293 0.1260
## district_citySubang -1.8303 0.1604 0.0827 0.3110
## district_citySukabumi -0.3954 0.6734 0.3613 1.2552
## district_citySukabumi_City -2.0110 0.1339 0.0684 0.2619
## district_citySumedang -1.7163 0.1797 0.0932 0.3468
## district_cityTasikmalaya -1.9532 0.1418 0.0727 0.2766
## district_cityTasikmalaya_City -1.7999 0.1653 0.0854 0.3201
## district_cityWest_Bandung -2.3210 0.0982 0.0493 0.1957Berdasarkan hasil model Negative Binomial, diperoleh beberapa interpretasi berikut:
Tren tahun: Nilai IRR sebesar 1,1990 menunjukkan bahwa tingkat laporan banjir meningkat sekitar 19,9% setiap tahun, dengan variabel lain dianggap tetap.
Bulan: Dibandingkan Januari sebagai bulan referensi, Februari memiliki tingkat laporan banjir 71,93% lebih tinggi. Sementara itu, Juni, Juli, Agustus, dan September masing-masing memiliki tingkat laporan lebih rendah sekitar 79,92%, 47,80%, 86,65%, dan 71,20%.
Wilayah: Kabupaten Bandung menjadi wilayah referensi. Sebagai contoh, IRR Kabupaten Bekasi sebesar 0,5224 menunjukkan tingkat laporan banjir sekitar 47,76% lebih rendah dibandingkan Kabupaten Bandung.
Prediksi Model
# Menghitung nilai prediksi
data_glm$prediksi <- predict(
model_nb,
type = "response"
)
# Melihat hasil aktual dan prediksi
head(data_glm[, c(
"district_city",
"tahun_tren",
"bulan",
"jumlah_banjir",
"prediksi"
)])## # A tibble: 6 × 5
## district_city tahun_tren bulan jumlah_banjir prediksi
## <fct> <dbl> <fct> <int> <dbl>
## 1 Bandung 0 1 3 5.42
## 2 Bandung 0 2 6 8.71
## 3 Bandung 0 3 3 6.94
## 4 Bandung 0 4 1 3.54
## 5 Bandung 0 5 1 3.95
## 6 Bandung 0 6 0 1.05Hasil menunjukkan perbandingan antara jumlah laporan banjir aktual dan hasil prediksi model Negative Binomial (NB) untuk Kabupaten Bandung pada tahun 2020.
Contohnya, pada Januari jumlah laporan aktual sebanyak 3, sedangkan model memprediksi sekitar 5,42 laporan. Pada Februari, jumlah aktual sebanyak 6, dengan prediksi 8,71 laporan.
Secara umum, nilai prediksi tidak harus sama persis dengan data aktual. Model digunakan untuk memperkirakan jumlah laporan banjir berdasarkan tren tahun, bulan, dan kabupaten/kota.
Visualisasi Aktual & Prediksi
ggplot(
data_glm,
aes(x = jumlah_banjir, y = prediksi)
) +
geom_point(alpha = 0.5) +
geom_abline(
intercept = 0,
slope = 1,
color = "red",
linetype = "dashed"
) +
labs(
title = "Perbandingan Jumlah Banjir Aktual dan Prediksi",
x = "Jumlah Laporan Banjir Aktual",
y = "Jumlah Laporan Banjir Prediksi"
) +
theme_minimal()
Grafik menunjukkan hubungan antara jumlah laporan banjir aktual (sumbu X) dan jumlah laporan hasil prediksi model NB (sumbu Y).
Titik hitam menunjukkan satu pengamatan, yaitu perbandingan antara jumlah laporan banjir aktual (sumbu X) dan jumlah laporan banjir hasil prediksi model NB (sumbu Y). Setiap titik mewakili satu data kabupaten/kota pada bulan dan tahun tertentu.
Garis merah putus-putus menunjukkan garis ideal, yaitu kondisi ketika nilai aktual sama dengan nilai prediksi (X = Y).
Cara membacanya: - Titik tepat pada garis merah: prediksi sama dengan nilai aktual. - Titik di atas garis merah: prediksi lebih tinggi daripada nilai aktual (overprediction). - Titik di bawah garis merah: prediksi lebih rendah daripada nilai aktual (underprediction). - Semakin dekat titik dengan garis merah, semakin baik ketepatan prediksinya.
Berdasarkan grafik, sebagian besar titik pengamatan belum berada dekat dengan garis merah putus-putus, yang menunjukkan bahwa nilai prediksi masih memiliki selisih dengan jumlah laporan banjir aktual. Pada beberapa pengamatan, nilai prediksi lebih tinggi daripada nilai aktual, sedangkan pada pengamatan lainnya nilai prediksi lebih rendah. Selain itu, terlihat beberapa titik dengan jumlah laporan aktual yang tinggi, tetapi nilai prediksinya relatif jauh lebih rendah.
Dengan demikian, model Negative Binomial belum mampu memprediksi jumlah laporan banjir secara akurat pada seluruh pengamatan. Meskipun demikian, model sudah dapat memberikan gambaran estimasi jumlah laporan banjir berdasarkan variabel yang digunakan.
KASUS 2. PEMODELAN JUMLAH LAPORAN KEJADIAN BANJIR MENGGUNAKAN REGRESI NONLINEAR
Selanjutnya menggunakan Generalized Additive Model (GAM) sebagai pendekatan regresi nonlinier untuk data jumlah laporan banjir.
# Memanggil paket
library(mgcv)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.# Membentuk model GAM Negative Binomial
model_gam <- gam(
jumlah_banjir ~ s(tahun_tren, k = 5) +
bulan + district_city +
offset(log(hari)),
family = nb(),
method = "REML",
data = data_glm
)
# Melihat hasil model
summary(model_gam)##
## Family: Negative Binomial(0.287)
## Link function: log
##
## Formula:
## jumlah_banjir ~ s(tahun_tren, k = 5) + bulan + district_city +
## offset(log(hari))
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2887 0.2757 -4.674 2.95e-06 ***
## bulan2 0.5421 0.2255 2.403 0.01624 *
## bulan3 0.2451 0.2268 1.081 0.27983
## bulan4 -0.3944 0.2328 -1.694 0.09027 .
## bulan5 -0.3184 0.2316 -1.375 0.16917
## bulan6 -1.6067 0.2564 -6.267 3.68e-10 ***
## bulan7 -0.6547 0.2357 -2.777 0.00548 **
## bulan8 -2.0222 0.2694 -7.506 6.10e-14 ***
## bulan9 -1.2476 0.2470 -5.051 4.40e-07 ***
## bulan10 -0.1491 0.2299 -0.648 0.51671
## bulan11 0.1085 0.2280 0.476 0.63424
## bulan12 0.1452 0.2275 0.638 0.52332
## district_cityBandung_City -2.2115 0.3549 -6.231 4.62e-10 ***
## district_cityBanjar_City -5.5437 0.7894 -7.023 2.17e-12 ***
## district_cityBekasi -0.6517 0.3268 -1.994 0.04613 *
## district_cityBekasi_City -1.0357 0.3309 -3.130 0.00175 **
## district_cityBogor 0.1698 0.3215 0.528 0.59745
## district_cityBogor_City -1.5046 0.3379 -4.452 8.50e-06 ***
## district_cityCiamis -1.2935 0.3344 -3.868 0.00011 ***
## district_cityCianjur -1.0539 0.3311 -3.183 0.00146 **
## district_cityCimahi_City -1.7749 0.3434 -5.169 2.35e-07 ***
## district_cityCirebon -0.6394 0.3267 -1.957 0.05034 .
## district_cityCirebon_City -1.9094 0.3465 -5.510 3.58e-08 ***
## district_cityDepok_City -1.5972 0.3397 -4.702 2.57e-06 ***
## district_cityGarut -0.7161 0.3274 -2.187 0.02873 *
## district_cityIndramayu -1.9620 0.3478 -5.641 1.69e-08 ***
## district_cityKarawang -0.2974 0.3240 -0.918 0.35878
## district_cityKuningan -2.3440 0.3592 -6.525 6.81e-11 ***
## district_cityMajalengka -1.9391 0.3472 -5.584 2.35e-08 ***
## district_cityPangandaran -1.7722 0.3433 -5.162 2.44e-07 ***
## district_cityPurwakarta -2.7982 0.3779 -7.404 1.32e-13 ***
## district_citySubang -1.8344 0.3447 -5.322 1.03e-07 ***
## district_citySukabumi -0.3926 0.3247 -1.209 0.22664
## district_citySukabumi_City -2.0133 0.3492 -5.766 8.12e-09 ***
## district_citySumedang -1.7179 0.3421 -5.022 5.13e-07 ***
## district_cityTasikmalaya -1.9403 0.3473 -5.587 2.31e-08 ***
## district_cityTasikmalaya_City -1.7965 0.3438 -5.225 1.74e-07 ***
## district_cityWest_Bandung -2.3231 0.3585 -6.480 9.20e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(tahun_tren) 1.006 1.013 39.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.126 Deviance explained = 29%
## -REML = 2908.4 Scale est. = 1 n = 1944Berdasarkan hasil pemodelan GAM Negative Binomial, diperoleh nilai deviance explained sebesar 29% dan adjusted R-squared sebesar 0,126. Artinya, model mampu menjelaskan sekitar 29% variasi jumlah laporan banjir melalui variabel yang digunakan.
Pada bagian smooth terms, variabel tahun_tren memiliki nilai p-value kurang dari 0,05, sehingga berpengaruh signifikan terhadap jumlah laporan banjir. Namun, nilai EDF sebesar 1,006 menunjukkan bahwa pola hubungan tahun dengan jumlah laporan banjir cenderung linear, bukan nonlinier yang kuat.
Sementara itu, berdasarkan koefisien parametrik, bulan Februari memiliki jumlah laporan banjir yang secara signifikan lebih tinggi dibandingkan Januari. Sebaliknya, bulan Juni, Juli, Agustus, dan September memiliki jumlah laporan yang lebih rendah dibandingkan Januari. Beberapa kabupaten/kota juga menunjukkan perbedaan signifikan dibandingkan Kabupaten Bandung sebagai kategori referensi.
Membandingkan GAM dan Negative Binomial
AIC(model_nb, model_gam)## df AIC
## model_nb 40.00000 5853.956
## model_gam 40.01255 5854.695Nilai AIC model Negative Binomial lebih kecil dibandingkan GAM, yaitu 5853,956. Selisih AIC kedua model sebesar 0,739, sehingga keduanya sebenarnya memiliki kualitas yang relatif sebanding. Namun, model Negative Binomial sedikit lebih unggul karena memiliki nilai AIC yang lebih rendah.
Kesimpulan: Model Negative Binomial dipilih sebagai model terbaik karena memiliki nilai AIC terendah. Selain itu, hasil GAM menunjukkan nilai EDF sebesar 1,006, yang mengindikasikan bahwa hubungan antara tahun dan jumlah laporan banjir cenderung linear. Dengan demikian, penggunaan regresi nonlinier belum memberikan peningkatan yang berarti dibandingkan model Negative Binomial sebelumnya.
KASUS 3. ANALISIS NON-LINEAR HUBUNGAN KONSENTRASI PM2.5 DAN PM10 MENGGUNAKAN REGRESI KUADRATIK PADA DATA KUALITAS UDARA WANLIU, BEIJING
Data dan Sumber Data
Data yang digunakan dalam analisis ini adalah data kualitas udara pada stasiun pemantauan Wanliu, Beijing, China, Chen, S. (2017). Beijing Multi-Site Air Quality [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C5RK5G.
Data tersebut merupakan data pengamatan kualitas udara setiap jam selama periode 1 Maret 2013 sampai 28 Februari 2017. Dataset berisi informasi mengenai konsentrasi polutan udara dan kondisi meteorologi.
Dalam analisis ini, variabel yang digunakan adalah:
PM2.5 Konsentrasi partikel udara berukuran ≤ 2,5 mikrometer PM10 Konsentrasi partikel udara berukuran ≤ 10 mikrometer
Dalam analisis ini, PM2.5 digunakan sebagai variabel respons (Y), sedangkan PM10 digunakan sebagai variabel prediktor (X). Analisis dilakukan untuk mengetahui bentuk hubungan antara konsentrasi PM2.5 dan PM10, khususnya kemungkinan adanya pola hubungan non-linear.
Persiapan Data
# Import data
data <- read.csv("C:/Users/rahay/Downloads/beijing+multi+site+air+quality+data/PRSA2017_Data_20130301-20170228/PRSA_Data_20130301-20170228/PRSA_Data_Wanliu_20130301-20170228.csv")
# Memilih variabel yang digunakan
data_nl <- data[, c("PM2.5", "PM10")]
# Menghapus data yang memiliki nilai kosong
data_nl <- na.omit(data_nl)
# Melihat jumlah data
nrow(data_nl)## [1] 34661sebanyak 34.661 observasi yang dapat digunakan dalam analisis.
Membentuk Model Non-Linear
Model regresi kuadratik dibentuk menggunakan sintaks berikut.
# Membentuk model regresi kuadratik
model_nl <- lm(PM2.5 ~ PM10 + I(PM10^2),
data = data_nl)
# Melihat hasil model
summary(model_nl)##
## Call:
## lm(formula = PM2.5 ~ PM10 + I(PM10^2), data = data_nl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -589.99 -15.24 3.66 16.22 265.02
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.065e+00 4.107e-01 -14.77 <2e-16 ***
## PM10 8.407e-01 5.132e-03 163.79 <2e-16 ***
## I(PM10^2) -1.665e-04 1.222e-05 -13.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.36 on 34658 degrees of freedom
## Multiple R-squared: 0.7799, Adjusted R-squared: 0.7799
## F-statistic: 6.14e+04 on 2 and 34658 DF, p-value: < 2.2e-16Berdasarkan hasil analisis, diperoleh persamaan regresi kuadratik sebagai berikut:
\[ \hat{Y} = -6.065 + 0.8407X - 0.0001665X^2 \]
Nilai koefisien determinasi yang diperoleh adalah (R^2 = 0.7799). Nilai tersebut menunjukkan bahwa sebesar 77,99% variasi konsentrasi PM2.5 dapat dijelaskan oleh PM10 melalui model regresi kuadratik, sedangkan 22,01% sisanya berkaitan dengan faktor lain yang tidak dimasukkan ke dalam model.
Koefisien PM10 yang bernilai positif menunjukkan bahwa pada bagian awal kurva, peningkatan PM10 berkaitan dengan peningkatan PM2.5. Sementara itu, koefisien PM10 kuadrat yang bernilai negatif menunjukkan bahwa pola hubungan berbentuk melengkung ke bawah. Dengan demikian, hubungan antara PM2.5 dan PM10 tidak hanya digambarkan sebagai hubungan linear.
Grafik Model
Visualisasi dilakukan untuk melihat pola hubungan antara PM2.5 dan PM10 serta membandingkannya dengan garis prediksi model regresi kuadratik.
# Scatter plot
plot(data_nl$PM10, data_nl$PM2.5,
main = "Analisis Non-Linear PM2.5 dan PM10",
xlab = "PM10",
ylab = "PM2.5",
pch = 16,
col = "gray")
# Membuat nilai prediksi
x <- seq(min(data_nl$PM10),
max(data_nl$PM10),
length.out = 500)
y <- predict(model_nl,
newdata = data.frame(PM10 = x))
# Menambahkan garis model
lines(x, y, col = "red", lwd = 2)
Berdasarkan grafik, konsentrasi PM2.5 cenderung meningkat seiring dengan meningkatnya konsentrasi PM10. Titik-titik berwarna abu-abu menunjukkan data pengamatan, sedangkan garis merah menunjukkan hasil prediksi model regresi kuadratik.
Garis model memperlihatkan pola yang sedikit melengkung ke bawah. Pola tersebut sejalan dengan koefisien PM10 kuadrat yang bernilai negatif. Meskipun demikian, sebaran titik pengamatan yang cukup luas menunjukkan bahwa konsentrasi PM2.5 tidak sepenuhnya dapat dijelaskan oleh PM10. Terdapat kemungkinan bahwa faktor lain, seperti kondisi meteorologi dan polutan lainnya, turut berkaitan dengan variasi konsentrasi PM2.5.
Menghitung Koefisien Determinasi
Koefisien determinasi digunakan untuk mengetahui proporsi variasi konsentrasi PM2.5 yang dapat dijelaskan oleh model regresi kuadratik.
# Menghitung R-squared
summary(model_nl)$r.squared## [1] 0.7798961diperoleh nilai \(R^2\) sebesar 0,7798961 atau 77,99%. Hal ini menunjukkan bahwa model regresi non-linear mampu menjelaskan 77,99% variasi konsentrasi PM2.5 berdasarkan PM10, sedangkan 22,01% sisanya dipengaruhi oleh faktor lain di luar model.