#===============================================================================
# 1. Library
#===============================================================================
# Menggabungkan semua library di awal agar bersih
libs <- c("gplots", "kableExtra", "tidyverse", "plm", "tinytex", "tseries",
"lmtest", "broom", "stargazer", "RColorBrewer", "corrplot", "car",
"readxl", "jsonlite", "GWmodel", "sf", "rnaturalearth", "patchwork",
"viridis")
lapply(libs, library, character.only = TRUE)
## Warning: package 'gplots' was built under R version 4.5.2
##
## ---------------------
## gplots 3.3.0 loaded:
## * Use citation('gplots') for citation info.
## * Homepage: https://talgalili.github.io/gplots/
## * Report issues: https://github.com/talgalili/gplots/issues
## * Ask questions: https://stackoverflow.com/questions/tagged/gplots
## * Suppress this message with: suppressPackageStartupMessages(library(gplots))
## ---------------------
##
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
##
## lowess
## Warning: package 'kableExtra' was built under R version 4.5.2
## Warning: package 'tidyverse' was built under R version 4.5.2
## Warning: package 'ggplot2' was built under R version 4.5.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Warning: package 'plm' was built under R version 4.5.2
##
## Attaching package: 'plm'
##
## The following objects are masked from 'package:dplyr':
##
## between, lag, lead
##
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Warning: package 'lmtest' was built under R version 4.5.2
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Warning: package 'stargazer' was built under R version 4.5.2
##
## Please cite as:
##
## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
##
## corrplot 0.95 loaded
## Warning: package 'car' was built under R version 4.5.2
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
##
##
## Attaching package: 'jsonlite'
##
## The following object is masked from 'package:purrr':
##
## flatten
## Warning: package 'GWmodel' was built under R version 4.5.2
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 4.5.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.5.2
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
## Warning: package 'sf' was built under R version 4.5.2
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
## Warning: package 'rnaturalearth' was built under R version 4.5.2
## Warning: package 'patchwork' was built under R version 4.5.2
## Warning: package 'viridis' was built under R version 4.5.2
## Loading required package: viridisLite
## [[1]]
## [1] "gplots" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "kableExtra" "gplots" "stats" "graphics" "grDevices"
## [6] "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [6] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [11] "kableExtra" "gplots" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "plm" "lubridate" "forcats" "stringr" "dplyr"
## [6] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [11] "tidyverse" "kableExtra" "gplots" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "tinytex" "plm" "lubridate" "forcats" "stringr"
## [6] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [11] "ggplot2" "tidyverse" "kableExtra" "gplots" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[6]]
## [1] "tseries" "tinytex" "plm" "lubridate" "forcats"
## [6] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [11] "tibble" "ggplot2" "tidyverse" "kableExtra" "gplots"
## [16] "stats" "graphics" "grDevices" "utils" "datasets"
## [21] "methods" "base"
##
## [[7]]
## [1] "lmtest" "zoo" "tseries" "tinytex" "plm"
## [6] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [11] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [16] "kableExtra" "gplots" "stats" "graphics" "grDevices"
## [21] "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "broom" "lmtest" "zoo" "tseries" "tinytex"
## [6] "plm" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "kableExtra" "gplots" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "stargazer" "broom" "lmtest" "zoo" "tseries"
## [6] "tinytex" "plm" "lubridate" "forcats" "stringr"
## [11] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [16] "ggplot2" "tidyverse" "kableExtra" "gplots" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"
##
## [[10]]
## [1] "RColorBrewer" "stargazer" "broom" "lmtest" "zoo"
## [6] "tseries" "tinytex" "plm" "lubridate" "forcats"
## [11] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [16] "tibble" "ggplot2" "tidyverse" "kableExtra" "gplots"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[11]]
## [1] "corrplot" "RColorBrewer" "stargazer" "broom" "lmtest"
## [6] "zoo" "tseries" "tinytex" "plm" "lubridate"
## [11] "forcats" "stringr" "dplyr" "purrr" "readr"
## [16] "tidyr" "tibble" "ggplot2" "tidyverse" "kableExtra"
## [21] "gplots" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[12]]
## [1] "car" "carData" "corrplot" "RColorBrewer" "stargazer"
## [6] "broom" "lmtest" "zoo" "tseries" "tinytex"
## [11] "plm" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "kableExtra" "gplots" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "readxl" "car" "carData" "corrplot" "RColorBrewer"
## [6] "stargazer" "broom" "lmtest" "zoo" "tseries"
## [11] "tinytex" "plm" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "kableExtra" "gplots" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[14]]
## [1] "jsonlite" "readxl" "car" "carData" "corrplot"
## [6] "RColorBrewer" "stargazer" "broom" "lmtest" "zoo"
## [11] "tseries" "tinytex" "plm" "lubridate" "forcats"
## [16] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [21] "tibble" "ggplot2" "tidyverse" "kableExtra" "gplots"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[15]]
## [1] "GWmodel" "Rcpp" "sp" "robustbase" "jsonlite"
## [6] "readxl" "car" "carData" "corrplot" "RColorBrewer"
## [11] "stargazer" "broom" "lmtest" "zoo" "tseries"
## [16] "tinytex" "plm" "lubridate" "forcats" "stringr"
## [21] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [26] "ggplot2" "tidyverse" "kableExtra" "gplots" "stats"
## [31] "graphics" "grDevices" "utils" "datasets" "methods"
## [36] "base"
##
## [[16]]
## [1] "sf" "GWmodel" "Rcpp" "sp" "robustbase"
## [6] "jsonlite" "readxl" "car" "carData" "corrplot"
## [11] "RColorBrewer" "stargazer" "broom" "lmtest" "zoo"
## [16] "tseries" "tinytex" "plm" "lubridate" "forcats"
## [21] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [26] "tibble" "ggplot2" "tidyverse" "kableExtra" "gplots"
## [31] "stats" "graphics" "grDevices" "utils" "datasets"
## [36] "methods" "base"
##
## [[17]]
## [1] "rnaturalearth" "sf" "GWmodel" "Rcpp"
## [5] "sp" "robustbase" "jsonlite" "readxl"
## [9] "car" "carData" "corrplot" "RColorBrewer"
## [13] "stargazer" "broom" "lmtest" "zoo"
## [17] "tseries" "tinytex" "plm" "lubridate"
## [21] "forcats" "stringr" "dplyr" "purrr"
## [25] "readr" "tidyr" "tibble" "ggplot2"
## [29] "tidyverse" "kableExtra" "gplots" "stats"
## [33] "graphics" "grDevices" "utils" "datasets"
## [37] "methods" "base"
##
## [[18]]
## [1] "patchwork" "rnaturalearth" "sf" "GWmodel"
## [5] "Rcpp" "sp" "robustbase" "jsonlite"
## [9] "readxl" "car" "carData" "corrplot"
## [13] "RColorBrewer" "stargazer" "broom" "lmtest"
## [17] "zoo" "tseries" "tinytex" "plm"
## [21] "lubridate" "forcats" "stringr" "dplyr"
## [25] "purrr" "readr" "tidyr" "tibble"
## [29] "ggplot2" "tidyverse" "kableExtra" "gplots"
## [33] "stats" "graphics" "grDevices" "utils"
## [37] "datasets" "methods" "base"
##
## [[19]]
## [1] "viridis" "viridisLite" "patchwork" "rnaturalearth"
## [5] "sf" "GWmodel" "Rcpp" "sp"
## [9] "robustbase" "jsonlite" "readxl" "car"
## [13] "carData" "corrplot" "RColorBrewer" "stargazer"
## [17] "broom" "lmtest" "zoo" "tseries"
## [21] "tinytex" "plm" "lubridate" "forcats"
## [25] "stringr" "dplyr" "purrr" "readr"
## [29] "tidyr" "tibble" "ggplot2" "tidyverse"
## [33] "kableExtra" "gplots" "stats" "graphics"
## [37] "grDevices" "utils" "datasets" "methods"
## [41] "base"
#===============================================================================
# 2. IMPORT DAN PERSIAPAN DATA
#===============================================================================
setwd("D:/Kuliah/Lomba/Law Phoria 2026")
banjir <- read_excel("data_deforestasi_final2.xlsx")
# Penamaan kolom agar konsisten dengan analisis selanjutnya
colnames(banjir) <- c("PROV", "TAHUN", "BNJR", "DEF", "KPDT", "CAC","SNGI")
banjir$PROV <- as.factor(toupper(banjir$PROV))
banjir$TAHUN <- as.factor(banjir$TAHUN)
# Cek Missing Value
print(colSums(is.na(banjir)))
## PROV TAHUN BNJR DEF KPDT CAC SNGI
## 0 0 0 0 0 0 0
#===============================================================================
# 3. EKSPLORASI DATA (VISUALISASI AWAL)
#===============================================================================
# Plot Korelasi antar variabel independen
# Penting untuk melihat hubungan awal sebelum regresi
# 1. Menghitung Matriks Korelasi
cor_matrix <- cor(banjir[, c("BNJR", "DEF", "KPDT", "CAC", "SNGI")], use = "complete.obs")
# 2. Visualisasi Corrplot dengan Estetika Custom
corrplot(
cor_matrix,
method = "color", # Menggunakan latar warna penuh pada kotak
type = "upper", # Hanya menampilkan bagian kanan atas
addCoef.col = "black", # Warna teks angka korelasi
number.cex = 0.75, # Ukuran teks angka korelasi
tl.cex = 0.9, # Ukuran teks label variabel (Top & Left)
tl.col = "gray15", # Warna label variabel
tl.srt = 0, # Sudut label (0 derajat agar mendatar)
tl.offset = 1,
col = colorRampPalette(c("#2166AC", "#FFFFFF", "#B2182B"))(200), # Gradien Biru-Putih-Merah
outline = TRUE, # Garis tepi pada setiap kotak
addgrid.col = "gray80", # Warna garis kisi-kisi
title = "Analisis Hubungan Antar Variabel (Correlation Matrix)",
mar = c(0, 0, 2, 0) # Memberi ruang untuk judul
)
# Visualisasi tren banjir per provinsi (Time Series Plot)
ggplot(data=banjir, aes(x=TAHUN, y=BNJR, group = PROV, colour = PROV)) +
geom_line(linewidth=1.2) +
geom_point(size=3, shape=19, fill="white") +
theme_bw() +
labs(title = "Dinamika Kejadian Banjir di Sumatera",
subtitle = "Periode 2018-2024",
x = "Tahun", y = "Jumlah Kejadian")
#===============================================================================
# 4. ANALISIS REGRESI DATA PANEL (GLOBAL)
#===============================================================================
# Uji Multikolinearitas (VIF)
check_vif <- lm(BNJR ~ DEF + KPDT + CAC + SNGI, data = banjir)
print("Nilai VIF (Harus < 10):")
## [1] "Nilai VIF (Harus < 10):"
print(vif(check_vif))
## DEF KPDT CAC SNGI
## 1.091874 1.073847 1.303694 1.290795
# 1. Common Effect Model (CEM)
cem <- plm(BNJR ~ DEF + KPDT + CAC + SNGI, data = banjir, model = "pooling")
# 2. Fixed Effect Model (FEM)
fem_ind <- plm(BNJR ~ DEF + KPDT + CAC + SNGI,
data = banjir, index = c("PROV", "TAHUN"),
model = "within", effect = "individual")
# 3. Fixed Effect Model (FEM) - TIME
fem_time <- plm(BNJR ~ DEF + KPDT + CAC + SNGI,
data = banjir, index = c("PROV", "TAHUN"),
model = "within", effect = "time")
# 4. Fixed Effect Model (REM) - TWO WAYS
fem_twoways <- plm(BNJR ~ DEF + KPDT + CAC + SNGI,
data = banjir, index = c("PROV", "TAHUN"),
model = "within", effect = "twoways")
#===============================================================================
# 5. TABEL PERBANDINGAN (LM TEST)
#===============================================================================
# 1. Jalankan Tes
test_ind <- plmtest(fem_twoways, type = "bp", effect = "individual")
test_time <- plmtest(fem_twoways, type = "bp", effect = "time")
test_two <- plmtest(fem_twoways, type = "bp", effect = "twoways")
# 2. Buat Data Frame Ringkas
tabel_lm <- data.frame(
Efek = c("Individual", "Time", "Twoways"),
BP_Stat = c(test_ind$statistic, test_time$statistic, test_two$statistic),
P_Value = c(test_ind$p.value, test_time$p.value, test_two$p.value)
)
# 3. Tambahkan Keputusan & Tampilkan Tabel
tabel_lm$Keputusan <- ifelse(tabel_lm$P_Value < 0.05, "Signifikan",
"Tidak Signifikan")
kable(tabel_lm, digits = 4, format = "markdown",
caption = "Hasil Uji LM Breusch-Pagan (Pemilihan Efek Model)")
Hasil Uji LM Breusch-Pagan (Pemilihan Efek Model)
| Individual |
0.3526 |
0.5526 |
Tidak Signifikan |
| Time |
2.3195 |
0.1278 |
Tidak Signifikan |
| Twoways |
2.6721 |
0.2629 |
Tidak Signifikan |
# 4. Hitung Perbandingan Performa Model (Ringkas)
mape <- function(actual, forecast) {
mean(abs((actual - forecast) / actual), na.rm = TRUE) * 100
}
lsdv_ind <- lm(BNJR ~ DEF + KPDT + CAC + SNGI + PROV, data = banjir)
lsdv_time <- lm(BNJR ~ DEF + KPDT + CAC + SNGI + TAHUN, data = banjir)
lsdv_twoways <- lm(BNJR ~ DEF + KPDT + CAC + SNGI + PROV + TAHUN, data = banjir)
compare <- data.frame(
Model = c("Individu", "Time", "Twoways"),
SSE = c(sum(fem_ind$residuals^2), sum(fem_time$residuals^2),
sum(fem_twoways$residuals^2)),
AIC = c(AIC(lsdv_ind), AIC(lsdv_time), AIC(lsdv_twoways)),
BIC = c(BIC(lsdv_ind), BIC(lsdv_time), BIC(lsdv_twoways)),
MAPE = c(mape(banjir$BNJR, predict(fem_ind)),
mape(banjir$BNJR, predict(fem_time)),
mape(banjir$BNJR, predict(fem_twoways))),
R2 = c(summary(fem_ind)$r.squared[1],
summary(fem_time)$r.squared[1],
summary(fem_twoways)$r.squared[1])
)
kable(compare, digits = 4, format = "markdown",
caption = "Perbandingan Kebaikan Model: Efek Individu, Time, dan Twoways")
Perbandingan Kebaikan Model: Efek Individu, Time, dan
Twoways
| Individu |
15075.20 |
604.7132 |
638.4406 |
277.6637 |
0.3136 |
| Time |
20131.48 |
618.9595 |
645.9415 |
291.3432 |
0.6578 |
| Twoways |
13887.84 |
610.9705 |
658.1889 |
285.0394 |
0.1043 |
# 5. Chow Test (CEM vs FEM Time)
chow_test <- pooltest(cem, fem_time)
tabel_chow <- data.frame(
Uji = "Chow Test (CEM vs FEM Time)",
Statistik = chow_test$statistic,
P_Value = chow_test$p.value,
Keputusan = ifelse(chow_test$p.value < 0.05, "Pilih FEM", "Pilih CEM")
)
kable(tabel_chow, digits = 5, format = "markdown",
caption = "Hasil Pemilihan Model: Chow Test")
Hasil Pemilihan Model: Chow Test
| F |
Chow Test (CEM vs FEM Time) |
2.29522 |
0.04658 |
Pilih FEM |
# 6. Tes BP pada Model REM
rem_gls <- plm(BNJR ~ DEF + KPDT + CAC + SNGI, data = banjir,
index = c("PROV", "TAHUN"),
effect = "time", model = "random", random.method = "nerlove")
rem_gls_ind <- plm(BNJR ~ DEF + KPDT + CAC + SNGI,
data = banjir,
index = c("PROV", "TAHUN"),
model = "random",
effect = "individual",
random.method = "nerlove")
bp_ind <- plmtest(rem_gls, type = "bp", effect = "individu")
bp_time <- plmtest(rem_gls, type = "bp", effect = "time")
bp_two <- plmtest(rem_gls, type = "bp", effect = "twoways")
tabel_rem_test <- data.frame(
Efek_Uji = c("Individu", "Time", "Twoways"),
BP_Statistik = c(bp_ind$statistic, bp_time$statistic, bp_two$statistic),
P_Value = c(bp_ind$p.value, bp_time$p.value, bp_two$p.value)
)
tabel_rem_test$Keputusan <- ifelse(tabel_rem_test$P_Value < 0.05,
"Signifikan",
"Tidak Signifikan")
kable(tabel_rem_test, digits = 4, format = "markdown",
caption = "Hasil Uji Lagrange Multiplier: Efek Random")
Hasil Uji Lagrange Multiplier: Efek Random
| Individu |
0.3526 |
0.5526 |
Tidak Signifikan |
| Time |
2.3195 |
0.1278 |
Tidak Signifikan |
| Twoways |
2.6721 |
0.2629 |
Tidak Signifikan |
# Hausman Test
h_test <- phtest(fem_time, rem_gls)
tabel_hausman <- data.frame(
Metode = "Hausman Test (FEM vs REM)",
Statistik = h_test$statistic,
P_Value = h_test$p.value,
Keputusan = ifelse(h_test$p.value < 0.05, "Pilih FEM", "Pilih REM")
)
kable(tabel_hausman, digits = 4, format = "markdown", caption = "Hasil Pemilihan Model: Hausman Test")
Hasil Pemilihan Model: Hausman Test
| chisq |
Hausman Test (FEM vs REM) |
0.8766 |
0.9279 |
Pilih REM |
# Uji Diagnostik
normality <- ks.test(rem_gls$residuals, "pnorm", mean(rem_gls$residuals), sd(rem_gls$residuals))
autocorr <- pbgtest(rem_gls)
homosced <- plmtest(rem_gls, type = "bp")
tabel_diagnostik <- data.frame(
Jenis_Uji = c("Normalitas (Kolmogorov-Smirnov)", "Autokorelasi (Breusch-Godfrey)", "Heteroskedastisitas (Breusch-Pagan)"),
Statistik = c(normality$statistic, autocorr$statistic, homosced$statistic),
P_Value = c(normality$p.value, autocorr$p.value, homosced$p.value)
)
tabel_diagnostik$Status_Asumsi <- ifelse(tabel_diagnostik$P_Value > 0.05,
"Terpenuhi",
"Terlanggar")
kable(tabel_diagnostik, digits = 4, format = "markdown", caption = "Ringkasan Uji Asumsi Klasik Model Panel")
Ringkasan Uji Asumsi Klasik Model Panel
| Normalitas (Kolmogorov-Smirnov) |
0.1108 |
0.3320 |
Terpenuhi |
| Autokorelasi (Breusch-Godfrey) |
12.8650 |
0.0755 |
Terpenuhi |
| Heteroskedastisitas (Breusch-Pagan) |
0.3526 |
0.5526 |
Terpenuhi |
#===============================================================================
# 6. PERSIAPAN DATA SPASIAL (INTEGRASI KOORDINAT)
#===============================================================================
# Mengambil data koordinat dari API GitHub (Yusuf Syaifudin)
# 1. Ambil data wilayah dari URL
wilayah_url <- "https://raw.githubusercontent.com/yusufsyaifudin/wilayah-indonesia/master/data/list_of_area/provinces.json"
wilayah_raw <- fromJSON(wilayah_url)
# 2. Standarisasi Nama Provinsi (Kapital)
banjir$PROV <- toupper(trimws(banjir$PROV))
wilayah_raw$name <- toupper(trimws(wilayah_raw$name))
# 3. Join data banjir dengan koordinat
data.sp.GWPR <- banjir %>%
left_join(wilayah_raw %>% select(name, latitude, longitude), by = c("PROV" = "name")) %>%
rename(LAT = latitude, LONG = longitude) %>%
# Pastikan koordinat adalah numerik (penting untuk gwr.basic)
mutate(LAT = as.numeric(LAT),
LONG = as.numeric(LONG)) %>%
# Mengatur ulang posisi kolom
select(PROV, TAHUN, LAT, LONG, everything())
coordinates(data.sp.GWPR) <- ~ LONG + LAT
#===============================================================================
# 7. PENENTUAN FUNGSI PEMBOBOT SPASIAL TERBAIK (GWPR)
#===============================================================================
# --- A. KERNEL ADAPTIVE (Bandwidth Berdasarkan Jumlah Tetangga) ---
# 1. Adaptive Bisquare
bwd.GWPR.bisquare.ad <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "bisquare", adaptive = TRUE)
## Adaptive bandwidth: 50 CV score: 29182.8
## Adaptive bandwidth: 39 CV score: 26316.59
## Adaptive bandwidth: 30 CV score: 27624.97
## Adaptive bandwidth: 42 CV score: 26316.59
hasil.GWPR.bisquare.ad <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.bisquare.ad, kernel = "bisquare", adaptive = TRUE)
# 2. Adaptive Gaussian
bwd.GWPR.gaussian.ad <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "gaussian", adaptive = TRUE)
## Adaptive bandwidth: 50 CV score: 28436.38
## Adaptive bandwidth: 39 CV score: 28406.31
## Adaptive bandwidth: 30 CV score: 28462.39
## Adaptive bandwidth: 42 CV score: 28406.31
hasil.GWPR.gaussian.ad <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.gaussian.ad, kernel = "gaussian", adaptive = TRUE)
# 3. Adaptive Exponential
bwd.GWPR.exponential.ad <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "exponential", adaptive = TRUE)
## Adaptive bandwidth: 50 CV score: 28176.23
## Adaptive bandwidth: 39 CV score: 27924.53
## Adaptive bandwidth: 30 CV score: 27717.7
## Adaptive bandwidth: 27 CV score: 27685.03
## Adaptive bandwidth: 22 CV score: 27685.03
hasil.GWPR.exponential.ad <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.exponential.ad, kernel = "exponential", adaptive = TRUE)
# --- B. KERNEL FIXED (Bandwidth Berdasarkan Jarak Tetap) ---
# 4. Fixed Bisquare
bwd.GWPR.bisquare.fix <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "bisquare", adaptive = FALSE)
## Fixed bandwidth: 7.873791 CV score: 29480.83
## Fixed bandwidth: 4.867243 CV score: 28871.67
## Fixed bandwidth: 3.009095 CV score: 293895.7
## Fixed bandwidth: 6.015642 CV score: 30431.17
## Fixed bandwidth: 4.157494 CV score: 29488.36
## Fixed bandwidth: 5.305893 CV score: 30126.75
## Fixed bandwidth: 4.596143 CV score: 28650.72
## Fixed bandwidth: 4.428594 CV score: 28862.68
## Fixed bandwidth: 4.699694 CV score: 28576.3
## Fixed bandwidth: 4.763692 CV score: 28631.32
## Fixed bandwidth: 4.660141 CV score: 28589.03
## Fixed bandwidth: 4.724139 CV score: 28586.82
## Fixed bandwidth: 4.684586 CV score: 28576.74
## Fixed bandwidth: 4.709031 CV score: 28578.71
## Fixed bandwidth: 4.693923 CV score: 28575.83
## Fixed bandwidth: 4.690357 CV score: 28575.94
## Fixed bandwidth: 4.696128 CV score: 28575.92
## Fixed bandwidth: 4.692561 CV score: 28575.84
## Fixed bandwidth: 4.694765 CV score: 28575.85
## Fixed bandwidth: 4.693403 CV score: 28575.83
## Fixed bandwidth: 4.693082 CV score: 28575.83
## Fixed bandwidth: 4.693602 CV score: 28575.83
## Fixed bandwidth: 4.69328 CV score: 28575.83
hasil.GWPR.bisquare.fix <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.bisquare.fix, kernel = "bisquare", adaptive = FALSE)
# 5. Fixed Gaussian
bwd.GWPR.gaussian.fix <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "gaussian", adaptive = FALSE)
## Fixed bandwidth: 7.873791 CV score: 28610.75
## Fixed bandwidth: 4.867243 CV score: 28537.84
## Fixed bandwidth: 3.009095 CV score: 28981.83
## Fixed bandwidth: 6.015642 CV score: 28531.73
## Fixed bandwidth: 6.725392 CV score: 28558.65
## Fixed bandwidth: 5.576993 CV score: 28522.82
## Fixed bandwidth: 5.305893 CV score: 28522.93
## Fixed bandwidth: 5.744542 CV score: 28525.15
## Fixed bandwidth: 5.473442 CV score: 28522.23
## Fixed bandwidth: 5.409444 CV score: 28522.24
## Fixed bandwidth: 5.512995 CV score: 28522.37
## Fixed bandwidth: 5.448997 CV score: 28522.2
## Fixed bandwidth: 5.433889 CV score: 28522.2
## Fixed bandwidth: 5.458334 CV score: 28522.2
## Fixed bandwidth: 5.443226 CV score: 28522.19
## Fixed bandwidth: 5.43966 CV score: 28522.19
## Fixed bandwidth: 5.44543 CV score: 28522.19
## Fixed bandwidth: 5.441864 CV score: 28522.19
hasil.GWPR.gaussian.fix <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.gaussian.fix, kernel = "gaussian", adaptive = FALSE)
# 6. Fixed Exponential
bwd.GWPR.exponential.fix <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
approach = "CV", kernel = "exponential", adaptive = FALSE)
## Fixed bandwidth: 7.873791 CV score: 28369
## Fixed bandwidth: 4.867243 CV score: 28158.94
## Fixed bandwidth: 3.009095 CV score: 27751.44
## Fixed bandwidth: 1.860696 CV score: 26731.4
## Fixed bandwidth: 1.150947 CV score: 26412.65
## Fixed bandwidth: 0.7122972 CV score: 38411.69
## Fixed bandwidth: 1.422047 CV score: 26090.55
## Fixed bandwidth: 1.589596 CV score: 26305.95
## Fixed bandwidth: 1.318496 CV score: 26051.53
## Fixed bandwidth: 1.254498 CV score: 26104.47
## Fixed bandwidth: 1.358049 CV score: 26052.28
## Fixed bandwidth: 1.294051 CV score: 26062.72
## Fixed bandwidth: 1.333604 CV score: 26049.31
## Fixed bandwidth: 1.342941 CV score: 26049.54
## Fixed bandwidth: 1.327833 CV score: 26049.76
## Fixed bandwidth: 1.33717 CV score: 26049.25
## Fixed bandwidth: 1.339374 CV score: 26049.31
## Fixed bandwidth: 1.335808 CV score: 26049.25
## Fixed bandwidth: 1.334966 CV score: 26049.27
## Fixed bandwidth: 1.336328 CV score: 26049.25
## Fixed bandwidth: 1.33665 CV score: 26049.25
## Fixed bandwidth: 1.336129 CV score: 26049.25
## Fixed bandwidth: 1.336451 CV score: 26049.25
hasil.GWPR.exponential.fix <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI, data = data.sp.GWPR,
bw = bwd.GWPR.exponential.fix, kernel = "exponential", adaptive = FALSE)
# 1. List Nama Model dan Objek Hasil untuk Iterasi (Opsional, tapi di sini kita buat manual agar scannable)
model_names <- c("Adaptive Bisquare", "Adaptive Gaussian", "Adaptive Exponential",
"Fixed Bisquare", "Fixed Gaussian", "Fixed Exponential")
# 2. Mengumpulkan Hasil Diagnostik dalam satu Data Frame
tabel_kernel <- data.frame(
Kernel = model_names,
RSS = c(hasil.GWPR.bisquare.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.gaussian.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.exponential.ad$GW.diagnostic$RSS.gw,
hasil.GWPR.bisquare.fix$GW.diagnostic$RSS.gw,
hasil.GWPR.gaussian.fix$GW.diagnostic$RSS.gw,
hasil.GWPR.exponential.fix$GW.diagnostic$RSS.gw),
R2 = c(hasil.GWPR.bisquare.ad$GW.diagnostic$gw.R2,
hasil.GWPR.gaussian.ad$GW.diagnostic$gw.R2,
hasil.GWPR.exponential.ad$GW.diagnostic$gw.R2,
hasil.GWPR.bisquare.fix$GW.diagnostic$gw.R2,
hasil.GWPR.gaussian.fix$GW.diagnostic$gw.R2,
hasil.GWPR.exponential.fix$GW.diagnostic$gw.R2),
R2Adj = c(hasil.GWPR.bisquare.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.gaussian.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.exponential.ad$GW.diagnostic$gwR2.adj,
hasil.GWPR.bisquare.fix$GW.diagnostic$gwR2.adj,
hasil.GWPR.gaussian.fix$GW.diagnostic$gwR2.adj,
hasil.GWPR.exponential.fix$GW.diagnostic$gwR2.adj),
AIC = c(hasil.GWPR.bisquare.ad$GW.diagnostic$AIC,
hasil.GWPR.gaussian.ad$GW.diagnostic$AIC,
hasil.GWPR.exponential.ad$GW.diagnostic$AIC,
hasil.GWPR.bisquare.fix$GW.diagnostic$AIC,
hasil.GWPR.gaussian.fix$GW.diagnostic$AIC,
hasil.GWPR.exponential.fix$GW.diagnostic$AIC),
AICc = c(hasil.GWPR.bisquare.ad$GW.diagnostic$AICc,
hasil.GWPR.gaussian.ad$GW.diagnostic$AICc,
hasil.GWPR.exponential.ad$GW.diagnostic$AICc,
hasil.GWPR.bisquare.fix$GW.diagnostic$AICc,
hasil.GWPR.gaussian.fix$GW.diagnostic$AICc,
hasil.GWPR.exponential.fix$GW.diagnostic$AICc)
)
min_aic <- min(tabel_kernel$AIC)
tabel_kernel$Keputusan <- ifelse(tabel_kernel$AIC == min_aic,
"Model Terbaik",
"-")
# 3. Tampilkan Tabel dengan kable
kable(tabel_kernel, digits = 3, format = "markdown",
caption = "Perbandingan Kebaikan Fungsi Pembobot Spasial (Kernel GWPR)")
Perbandingan Kebaikan Fungsi Pembobot Spasial (Kernel
GWPR)
| Adaptive Bisquare |
17775.72 |
0.728 |
0.640 |
601.282 |
628.630 |
- |
| Adaptive Gaussian |
22440.74 |
0.656 |
0.609 |
609.656 |
621.168 |
- |
| Adaptive Exponential |
20412.53 |
0.687 |
0.622 |
605.660 |
621.709 |
- |
| Fixed Bisquare |
17978.67 |
0.725 |
0.636 |
602.146 |
629.664 |
- |
| Fixed Gaussian |
22982.44 |
0.648 |
0.605 |
610.786 |
621.445 |
- |
| Fixed Exponential |
17094.97 |
0.738 |
0.647 |
598.705 |
626.437 |
Model Terbaik |
#===============================================================================
# 6. MODELING GWPR (LOOPING PER TAHUN) - PERBAIKAN
#===============================================================================
list_sdf_tahunan <- list()
tahun_analisis <- sort(unique(data.sp.GWPR$TAHUN))
for (th in tahun_analisis) {
# 1 & 2. Filter data per tahun DAN pastikan formatnya SpatialPointsDataFrame
# Kita ambil dari data asli sebelum konversi jika ada, atau gunakan index:
data_tahun_ini <- data.sp.GWPR[data.sp.GWPR$TAHUN == th, ]
# 3. Optimasi Bandwidth Fixed Bisquare (karena adaptive = FALSE)
bw_optimum <- bw.gwr(BNJR ~ DEF + KPDT + CAC + SNGI,
data = data_tahun_ini,
approach = "CV",
kernel = "exponential",
adaptive = FALSE)
# 4. Estimasi Model GWR
model_gwr <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI,
data = data_tahun_ini,
bw = bw_optimum,
kernel = "exponential",
adaptive = FALSE)
# 5. Ekstraksi hasil SDF ke Data Frame
hasil_sdf <- as.data.frame(model_gwr$SDF)
hasil_sdf$TAHUN <- th
hasil_sdf$PROV <- data_tahun_ini$PROV
# Masukkan ke list
list_sdf_tahunan[[as.character(th)]] <- hasil_sdf
message(paste("Selesai memproses tahun:", th))
}
## Fixed bandwidth: 7.873791 CV score: 2078.135
## Fixed bandwidth: 4.867243 CV score: 1863.494
## Fixed bandwidth: 3.009095 CV score: 1587.763
## Fixed bandwidth: 1.860696 CV score: 1468.506
## Fixed bandwidth: 1.150947 CV score: 1998.253
## Fixed bandwidth: 2.299345 CV score: 1472.655
## Fixed bandwidth: 1.589596 CV score: 1549.847
## Fixed bandwidth: 2.028245 CV score: 1456.944
## Fixed bandwidth: 2.131796 CV score: 1459.079
## Fixed bandwidth: 1.964247 CV score: 1458.843
## Fixed bandwidth: 2.067798 CV score: 1457.068
## Fixed bandwidth: 2.0038 CV score: 1457.346
## Fixed bandwidth: 2.043353 CV score: 1456.882
## Fixed bandwidth: 2.05269 CV score: 1456.913
## Fixed bandwidth: 2.037582 CV score: 1456.89
## Fixed bandwidth: 2.04692 CV score: 1456.888
## Fixed bandwidth: 2.041149 CV score: 1456.883
## Fixed bandwidth: 2.044715 CV score: 1456.884
## Fixed bandwidth: 2.042511 CV score: 1456.882
## Fixed bandwidth: 2.041991 CV score: 1456.882
## Selesai memproses tahun: 2018
## Fixed bandwidth: 7.873791 CV score: 1843.687
## Fixed bandwidth: 4.867243 CV score: 1978.417
## Fixed bandwidth: 9.731939 CV score: 1806.714
## Fixed bandwidth: 10.88034 CV score: 1790.894
## Fixed bandwidth: 11.59009 CV score: 1782.855
## Fixed bandwidth: 12.02874 CV score: 1778.412
## Fixed bandwidth: 12.29984 CV score: 1775.841
## Fixed bandwidth: 12.46739 CV score: 1774.313
## Fixed bandwidth: 12.57094 CV score: 1773.392
## Fixed bandwidth: 12.63494 CV score: 1772.831
## Fixed bandwidth: 12.67449 CV score: 1772.487
## Fixed bandwidth: 12.69893 CV score: 1772.276
## Fixed bandwidth: 12.71404 CV score: 1772.146
## Fixed bandwidth: 12.72338 CV score: 1772.065
## Fixed bandwidth: 12.72915 CV score: 1772.016
## Fixed bandwidth: 12.73272 CV score: 1771.985
## Fixed bandwidth: 12.73492 CV score: 1771.966
## Fixed bandwidth: 12.73628 CV score: 1771.955
## Fixed bandwidth: 12.73712 CV score: 1771.947
## Fixed bandwidth: 12.73764 CV score: 1771.943
## Fixed bandwidth: 12.73797 CV score: 1771.94
## Fixed bandwidth: 12.73816 CV score: 1771.938
## Fixed bandwidth: 12.73829 CV score: 1771.937
## Fixed bandwidth: 12.73836 CV score: 1771.937
## Fixed bandwidth: 12.73841 CV score: 1771.936
## Fixed bandwidth: 12.73844 CV score: 1771.936
## Selesai memproses tahun: 2019
## Fixed bandwidth: 7.873791 CV score: 3615.306
## Fixed bandwidth: 4.867243 CV score: 3813.625
## Fixed bandwidth: 9.731939 CV score: 3573.574
## Fixed bandwidth: 10.88034 CV score: 3557.843
## Fixed bandwidth: 11.59009 CV score: 3550.389
## Fixed bandwidth: 12.02874 CV score: 3546.433
## Fixed bandwidth: 12.29984 CV score: 3544.199
## Fixed bandwidth: 12.46739 CV score: 3542.891
## Fixed bandwidth: 12.57094 CV score: 3542.109
## Fixed bandwidth: 12.63494 CV score: 3541.635
## Fixed bandwidth: 12.67449 CV score: 3541.346
## Fixed bandwidth: 12.69893 CV score: 3541.169
## Fixed bandwidth: 12.71404 CV score: 3541.059
## Fixed bandwidth: 12.72338 CV score: 3540.992
## Fixed bandwidth: 12.72915 CV score: 3540.951
## Fixed bandwidth: 12.73272 CV score: 3540.925
## Fixed bandwidth: 12.73492 CV score: 3540.909
## Fixed bandwidth: 12.73628 CV score: 3540.9
## Fixed bandwidth: 12.73712 CV score: 3540.894
## Fixed bandwidth: 12.73764 CV score: 3540.89
## Fixed bandwidth: 12.73797 CV score: 3540.888
## Fixed bandwidth: 12.73816 CV score: 3540.886
## Fixed bandwidth: 12.73829 CV score: 3540.885
## Fixed bandwidth: 12.73836 CV score: 3540.885
## Fixed bandwidth: 12.73841 CV score: 3540.884
## Fixed bandwidth: 12.73844 CV score: 3540.884
## Selesai memproses tahun: 2020
## Fixed bandwidth: 7.873791 CV score: 8570.033
## Fixed bandwidth: 4.867243 CV score: 9158.03
## Fixed bandwidth: 9.731939 CV score: 8470.773
## Fixed bandwidth: 10.88034 CV score: 8436.059
## Fixed bandwidth: 11.59009 CV score: 8420.176
## Fixed bandwidth: 12.02874 CV score: 8411.902
## Fixed bandwidth: 12.29984 CV score: 8407.277
## Fixed bandwidth: 12.46739 CV score: 8404.587
## Fixed bandwidth: 12.57094 CV score: 8402.984
## Fixed bandwidth: 12.63494 CV score: 8402.015
## Fixed bandwidth: 12.67449 CV score: 8401.425
## Fixed bandwidth: 12.69893 CV score: 8401.063
## Fixed bandwidth: 12.71404 CV score: 8400.841
## Fixed bandwidth: 12.72338 CV score: 8400.703
## Fixed bandwidth: 12.72915 CV score: 8400.619
## Fixed bandwidth: 12.73272 CV score: 8400.567
## Fixed bandwidth: 12.73492 CV score: 8400.535
## Fixed bandwidth: 12.73628 CV score: 8400.515
## Fixed bandwidth: 12.73712 CV score: 8400.502
## Fixed bandwidth: 12.73764 CV score: 8400.495
## Fixed bandwidth: 12.73797 CV score: 8400.49
## Fixed bandwidth: 12.73816 CV score: 8400.487
## Fixed bandwidth: 12.73829 CV score: 8400.485
## Fixed bandwidth: 12.73836 CV score: 8400.484
## Fixed bandwidth: 12.73841 CV score: 8400.484
## Fixed bandwidth: 12.73844 CV score: 8400.483
## Selesai memproses tahun: 2021
## Fixed bandwidth: 7.873791 CV score: 8740.078
## Fixed bandwidth: 4.867243 CV score: 9775.692
## Fixed bandwidth: 9.731939 CV score: 8514.005
## Fixed bandwidth: 10.88034 CV score: 8426.11
## Fixed bandwidth: 11.59009 CV score: 8383.592
## Fixed bandwidth: 12.02874 CV score: 8360.726
## Fixed bandwidth: 12.29984 CV score: 8347.705
## Fixed bandwidth: 12.46739 CV score: 8340.043
## Fixed bandwidth: 12.57094 CV score: 8335.446
## Fixed bandwidth: 12.63494 CV score: 8332.657
## Fixed bandwidth: 12.67449 CV score: 8330.952
## Fixed bandwidth: 12.69893 CV score: 8329.906
## Fixed bandwidth: 12.71404 CV score: 8329.262
## Fixed bandwidth: 12.72338 CV score: 8328.865
## Fixed bandwidth: 12.72915 CV score: 8328.62
## Fixed bandwidth: 12.73272 CV score: 8328.469
## Fixed bandwidth: 12.73492 CV score: 8328.376
## Fixed bandwidth: 12.73628 CV score: 8328.318
## Fixed bandwidth: 12.73712 CV score: 8328.282
## Fixed bandwidth: 12.73764 CV score: 8328.26
## Fixed bandwidth: 12.73797 CV score: 8328.247
## Fixed bandwidth: 12.73816 CV score: 8328.238
## Fixed bandwidth: 12.73829 CV score: 8328.233
## Fixed bandwidth: 12.73836 CV score: 8328.23
## Fixed bandwidth: 12.73841 CV score: 8328.228
## Fixed bandwidth: 12.73844 CV score: 8328.227
## Selesai memproses tahun: 2022
## Fixed bandwidth: 7.873791 CV score: 2427.261
## Fixed bandwidth: 4.867243 CV score: 2762.414
## Fixed bandwidth: 9.731939 CV score: 2338.868
## Fixed bandwidth: 10.88034 CV score: 2301.82
## Fixed bandwidth: 11.59009 CV score: 2283.21
## Fixed bandwidth: 12.02874 CV score: 2272.991
## Fixed bandwidth: 12.29984 CV score: 2267.101
## Fixed bandwidth: 12.46739 CV score: 2263.61
## Fixed bandwidth: 12.57094 CV score: 2261.507
## Fixed bandwidth: 12.63494 CV score: 2260.227
## Fixed bandwidth: 12.67449 CV score: 2259.444
## Fixed bandwidth: 12.69893 CV score: 2258.963
## Fixed bandwidth: 12.71404 CV score: 2258.666
## Fixed bandwidth: 12.72338 CV score: 2258.483
## Fixed bandwidth: 12.72915 CV score: 2258.371
## Fixed bandwidth: 12.73272 CV score: 2258.301
## Fixed bandwidth: 12.73492 CV score: 2258.258
## Fixed bandwidth: 12.73628 CV score: 2258.231
## Fixed bandwidth: 12.73712 CV score: 2258.215
## Fixed bandwidth: 12.73764 CV score: 2258.205
## Fixed bandwidth: 12.73797 CV score: 2258.199
## Fixed bandwidth: 12.73816 CV score: 2258.195
## Fixed bandwidth: 12.73829 CV score: 2258.192
## Fixed bandwidth: 12.73836 CV score: 2258.191
## Fixed bandwidth: 12.73841 CV score: 2258.19
## Fixed bandwidth: 12.73844 CV score: 2258.189
## Selesai memproses tahun: 2023
## Fixed bandwidth: 7.873791 CV score: 9403.702
## Fixed bandwidth: 4.867243 CV score: 7857.661
## Fixed bandwidth: 3.009095 CV score: 6210.329
## Fixed bandwidth: 1.860696 CV score: 7465.143
## Fixed bandwidth: 3.718844 CV score: 6886.084
## Fixed bandwidth: 2.570446 CV score: 5924.863
## Fixed bandwidth: 2.299345 CV score: 5982.513
## Fixed bandwidth: 2.737995 CV score: 6002.29
## Fixed bandwidth: 2.466895 CV score: 5911.925
## Fixed bandwidth: 2.402896 CV score: 5923
## Fixed bandwidth: 2.506448 CV score: 5912.863
## Fixed bandwidth: 2.442449 CV score: 5914.187
## Fixed bandwidth: 2.482002 CV score: 5911.641
## Fixed bandwidth: 2.49134 CV score: 5911.869
## Fixed bandwidth: 2.476232 CV score: 5911.652
## Fixed bandwidth: 2.485569 CV score: 5911.692
## Fixed bandwidth: 2.479798 CV score: 5911.631
## Fixed bandwidth: 2.478436 CV score: 5911.634
## Fixed bandwidth: 2.48064 CV score: 5911.633
## Fixed bandwidth: 2.479278 CV score: 5911.632
## Fixed bandwidth: 2.48012 CV score: 5911.632
## Fixed bandwidth: 2.479599 CV score: 5911.631
## Selesai memproses tahun: 2024
# Gabungkan semua hasil menjadi satu dataframe besar untuk langkah prediksi
gwpr_results_full <- do.call(rbind, list_sdf_tahunan)
#===============================================================================
# 7. ANALISIS SIGNIFIKANSI PARAMETER LOKAL
#===============================================================================
# Menentukan tingkat kepercayaan 95% (Alpha 0.05)
# Nilai t-kritis standar adalah 1.96
t_kritis_val <- 1.96
gwpr_summary_sig <- gwpr_results_full %>%
mutate(
DEF_sig = ifelse(abs(DEF_TV) >= t_kritis_val, "Signifikan", "Tidak"),
KPDT_sig = ifelse(abs(KPDT_TV) >= t_kritis_val, "Signifikan", "Tidak"),
CAC_sig = ifelse(abs(CAC_TV) >= t_kritis_val, "Signifikan", "Tidak"),
SNGI_sig = ifelse(abs(SNGI_TV) >= t_kritis_val, "Signifikan", "Tidak")
)
#===============================================================================
# 8. OUTPUT TABEL PEMBAHASAN
#===============================================================================
# Tabel A: Ringkasan Signifikansi Per Tahun
tabel_tahun_sig <- gwpr_summary_sig %>%
group_by(TAHUN) %>%
summarise(
Deforestasi = sum(DEF_sig == "Signifikan"),
Kepadatan_Pend = sum(KPDT_sig == "Signifikan"),
Curah_Hujan = sum(CAC_sig == "Signifikan"),
Jarak_Sungai = sum(SNGI_sig == "Signifikan")
)
kable(tabel_tahun_sig, caption = "Tabel 1. Perkembangan Jumlah Provinsi Signifikan (2018-2024)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Tabel 1. Perkembangan Jumlah Provinsi Signifikan (2018-2024)
|
TAHUN
|
Deforestasi
|
Kepadatan_Pend
|
Curah_Hujan
|
Jarak_Sungai
|
|
2018
|
0
|
0
|
10
|
4
|
|
2019
|
0
|
8
|
10
|
0
|
|
2020
|
0
|
0
|
10
|
10
|
|
2021
|
0
|
0
|
8
|
10
|
|
2022
|
0
|
0
|
2
|
10
|
|
2023
|
10
|
0
|
10
|
9
|
|
2024
|
0
|
0
|
0
|
0
|
# Tabel B: Daftar Provinsi Signifikan Positif/Negatif
tahun_list <- sort(unique(gwpr_summary_sig$TAHUN))
vars_to_check <- c("DEF", "KPDT", "CAC", "SNGI")
list_rekap_tahunan <- list()
for (th in tahun_list) {
data_tahunan <- gwpr_summary_sig %>% filter(TAHUN == th)
rekap_temp <- data.frame(
Tahun = th,
Arah = c("Positif Signifikan", "Negatif Signifikan")
)
for (v in vars_to_check) {
pos <- data_tahunan$PROV[data_tahunan[[paste0(v, "_sig")]] == "Signifikan" & data_tahunan[[v]] > 0]
neg <- data_tahunan$PROV[data_tahunan[[paste0(v, "_sig")]] == "Signifikan" & data_tahunan[[v]] < 0]
rekap_temp[[v]] <- c(
if(length(pos) > 0) paste(sort(pos), collapse = ", ") else "-",
if(length(neg) > 0) paste(sort(neg), collapse = ", ") else "-"
)
}
list_rekap_tahunan[[as.character(th)]] <- rekap_temp
}
tabel_rekap_final <- do.call(rbind, list_rekap_tahunan)
kable(tabel_rekap_final,
format = "markdown",
row.names = FALSE,
caption = "Tabel 2. Evolusi Klasifikasi Pengaruh Variabel per Provinsi (2018-2024)")
Tabel 2. Evolusi Klasifikasi Pengaruh Variabel per Provinsi
(2018-2024)
| 2018 |
Positif Signifikan |
- |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
- |
| 2018 |
Negatif Signifikan |
- |
- |
- |
ACEH, RIAU, SUMATERA BARAT, SUMATERA UTARA |
| 2019 |
Positif Signifikan |
- |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
- |
| 2019 |
Negatif Signifikan |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN RIAU, RIAU, SUMATERA
BARAT, SUMATERA SELATAN, SUMATERA UTARA |
- |
- |
| 2020 |
Positif Signifikan |
- |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
| 2020 |
Negatif Signifikan |
- |
- |
- |
- |
| 2021 |
Positif Signifikan |
- |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, SUMATERA SELATAN, SUMATERA UTARA |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
| 2021 |
Negatif Signifikan |
- |
- |
- |
- |
| 2022 |
Positif Signifikan |
- |
- |
ACEH, SUMATERA UTARA |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
| 2022 |
Negatif Signifikan |
- |
- |
- |
- |
| 2023 |
Positif Signifikan |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA SELATAN, SUMATERA UTARA |
| 2023 |
Negatif Signifikan |
- |
- |
- |
- |
| 2024 |
Positif Signifikan |
- |
- |
- |
- |
| 2024 |
Negatif Signifikan |
- |
- |
- |
- |
#===============================================================================
# 9. VISUALISASI PETA SPASIAL FINAL
#===============================================================================
data_2024_akhir <- gwpr_summary_sig %>% filter(TAHUN == 2024)
# Load peta Sumatera
indo_map <- ne_states(country = "indonesia", returnclass = "sf")
sumatera_map <- indo_map %>% filter(name %in% c("Aceh", "Sumatera Utara", "Sumatera Barat", "Riau", "Jambi",
"Sumatera Selatan", "Bengkulu", "Lampung",
"Kepulauan Bangka Belitung", "Kepulauan Riau"))
# Convert hasil akhir ke SF object untuk pemetaan
peta_sf_2024 <- st_as_sf(data_2024_akhir, coords = c("LONG", "LAT"), crs = 4326)
# Peta Sebaran Koefisien DEF dengan Background Peta Asli
ggplot() +
geom_sf(data = sumatera_map, fill = "gray95", color = "gray80") +
geom_sf(data = peta_sf_2024, aes(color = DEF, size = Local_R2), alpha = 0.8) +
scale_color_viridis_c(option = "plasma", name = "Koefisien DEF") +
scale_size_continuous(range = c(8, 16), name = "Akurasi (R2)") +
geom_sf_text(data = peta_sf_2024, aes(label = PROV), size = 3, vjust = 4, fontface = "bold") +
labs(title = "Peta Sebaran Dampak Deforestasi terhadap Banjir",
subtitle = "Analisis GWPR Tahun 2024",
caption = "Warna menunjukkan kekuatan pengaruh | Ukuran menunjukkan akurasi lokal") +
theme_minimal()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# 1. Pastikan data sudah menggabungkan koordinat dan identitas tahun
# Kita gunakan data parameter.GWPR yang sudah memiliki 70 baris (10 prov x 7 thn)
gwpr_all_years_sf <- st_as_sf(gwpr_summary_sig, coords = c("LONG", "LAT"), crs = 4326)
# 1. Hitung nilai min dan max dari koefisien DEF untuk label legenda
min_val <- min(gwpr_summary_sig$DEF, na.rm = TRUE)
max_val <- max(gwpr_summary_sig$DEF, na.rm = TRUE)
# 2. Plotting
ggplot() +
geom_sf(data = sumatera_map, fill = "gray97", color = "gray80", size = 0.2) +
# PERBAIKAN: Pastikan color diarahkan ke DEF jika ingin melihat dampak deforestasi
geom_sf(data = gwpr_all_years_sf, aes(color = DEF, size = Local_R2), alpha = 0.8) +
# PERBAIKAN: Gunakan nilai breaks yang sesuai dengan data DEF
scale_color_viridis_c(option = "plasma",
name = "Dampak Deforestasi",
breaks = c(min_val, max_val),
labels = c("Terendah", "Tertinggi")) +
scale_size_continuous(range = c(4, 8), name = "Local R2") +
facet_wrap(~TAHUN, ncol = 4) +
labs(title = "Evolusi Spasial Dampak Deforestasi terhadap Banjir di Sumatera",
subtitle = "Analisis Geographically Weighted Panel Regression (2018-2024)",
caption = "Warna Titik: Koefisien DEF | Ukuran Titik: Local R2") +
theme_minimal() +
theme(
strip.background = element_rect(fill = "gray20"),
strip.text = element_text(color = "white", face = "bold"),
legend.position = "bottom", # Pastikan posisi legenda benar
axis.text = element_blank(),
axis.ticks = element_blank()
)
#===============================================================================
# 10. PREDIKSI DINAMIS (MEMANFAATKAN DATA 2018-2024)
#===============================================================================
# A. Membuat data rangka untuk tahun 2025-2030
tahun_proyeksi <- 2025:2030
provinsi_list <- unique(banjir$PROV)
data_proyeksi <- expand.grid(PROV = provinsi_list, TAHUN = tahun_proyeksi)
# B. Prediksi nilai Variabel Independen (DEF, KPDT, dll) menggunakan Tren Linear
# Kita asumsikan tren deforestasi dan variabel lain berlanjut secara linear
predict_trend <- function(df, target_prov, target_var) {
sub_data <- banjir %>% filter(PROV == target_prov)
sub_data$TAHUN_NUM <- as.numeric(as.character(sub_data$TAHUN))
model_trend <- lm(as.formula(paste(target_var, "~ TAHUN_NUM")), data = sub_data)
return(predict(model_trend, newdata = data.frame(TAHUN_NUM = 2025:2030)))
}
# Mengisi nilai prediksi variabel independen ke data_proyeksi
for(p in provinsi_list) {
for(v in c("DEF", "KPDT", "CAC", "SNGI")) {
data_proyeksi[data_proyeksi$PROV == p, v] <- predict_trend(banjir, p, v)
}
}
# A. Forecasting Koefisien GWPR per Provinsi (Intercept, DEF, KPDT, dll)
# Kita hitung tren perubahan pengaruh variabel dari tahun ke tahun
predict_coef_trend <- function(target_prov, target_coef) {
sub_coef <- gwpr_results_full %>% filter(PROV == target_prov)
sub_coef$TAHUN_NUM <- as.numeric(as.character(sub_coef$TAHUN))
# Model tren untuk koefisien itu sendiri
model_c <- lm(as.formula(paste(target_coef, "~ TAHUN_NUM")), data = sub_coef)
return(predict(model_c, newdata = data.frame(TAHUN_NUM = 2025:2030)))
}
# B. Buat dataframe proyeksi baru
tahun_fut <- 2025:2030
prov_list <- unique(banjir$PROV)
data_future_dynamic <- expand.grid(PROV = prov_list, TAHUN = tahun_fut)
# C. Isi Koefisien Prediksi & Variabel Prediksi (Looping Dinamis)
for(p in prov_list) {
# Prediksi Koefisien (Intercept + Variabel)
for(co in c("Intercept", "DEF", "KPDT", "CAC", "SNGI")) {
data_future_dynamic[data_future_dynamic$PROV == p, paste0(co, "_pred")] <- predict_coef_trend(p, co)
}
# Prediksi Nilai Variabel Independen (X)
for(v in c("DEF", "KPDT", "CAC", "SNGI")) {
data_future_dynamic[data_future_dynamic$PROV == p, v] <- predict_trend(banjir, p, v)
}
}
# D. Hitung BNJR dengan Koefisien yang Berubah Tiap Tahun
# Ditambah sedikit noise stokastik agar visualisasi memiliki variabilitas alami
set.seed(123)
data_future_dynamic <- data_future_dynamic %>%
mutate(
BNJR = Intercept_pred +
(DEF * DEF_pred) +
(KPDT * KPDT_pred) +
(CAC * CAC_pred) +
(SNGI * SNGI_pred),
# Menambahkan variabilitas acak (noise) berdasarkan standar deviasi historis
Noise = rnorm(n(), mean = 0, sd = sd(banjir$BNJR) * 0.2),
BNJR = pmax(0, BNJR + Noise) # Pastikan tidak negatif
)
#===============================================================================
# 11. VISUALISASI TIME SERIES REVISI (NAIK-TURUN)
#===============================================================================
# 1. Pastikan data sudah siap
data_viz <- rbind(
banjir %>% select(PROV, TAHUN, BNJR) %>%
mutate(Tipe = "Historis", TAHUN = as.numeric(as.character(TAHUN))),
data_future_dynamic %>% select(PROV, TAHUN, BNJR) %>%
mutate(Tipe = "Proyeksi GWPR", TAHUN = as.numeric(as.character(TAHUN)))
)
# 2. Buat Plot dengan membagi geom_line
p_timeseries <- ggplot(data_viz, aes(x = TAHUN, y = BNJR, color = PROV, group = PROV)) +
# Area Prediksi (Background)
annotate("rect", xmin = 2024.5, xmax = 2030.5, ymin = -Inf, ymax = Inf,
fill = "gray90", alpha = 0.4) +
# LAYER 1: Garis Historis (2018-2024)
geom_line(data = subset(data_viz, Tipe == "Historis"),
linewidth = 1, linetype = "solid") +
# LAYER 2: Garis Proyeksi (2024-2030)
# Kita tambahkan tahun 2024 ke data proyeksi agar garisnya menyambung sempurna
geom_line(data = subset(data_viz, Tipe == "Proyeksi GWPR" | TAHUN == 2024),
linewidth = 1) +
# LAYER 3: Titik-titik data
geom_point(aes(shape = Tipe), size = 2, alpha = 0.7) +
# Garis Vertikal Pemisah
geom_vline(xintercept = 2024.5, linetype = "dotted", color = "black", alpha = 0.5) +
# Pengaturan Skala dan Label
scale_x_continuous(breaks = 2018:2030) +
labs(title = "Evolusi & Proyeksi Dinamis Kejadian Banjir di Sumatera",
subtitle = "Estimasi Dinamis: Tren Koefisien Lokal GWPR 2018-2030",
x = "Tahun", y = "Jumlah Kejadian Banjir") +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text.x = element_text(angle = 45, vjust = 0.5),
panel.grid.minor = element_blank()
)
# Tampilkan Plot
print(p_timeseries)
#===============================================================================
# 12. SKENARIO KEBIJAKAN (2025-2030)
#===============================================================================
# Fungsi dasar untuk menghitung BNJR berdasarkan skenario
calc_flood <- function(df, noise_level = 0.15) {
set.seed(42)
df %>% mutate(
BNJR = Intercept_pred + (DEF * DEF_pred) + (KPDT * KPDT_pred) +
(CAC * CAC_pred) + (SNGI * SNGI_pred),
Noise = rnorm(n(), mean = 0, sd = sd(banjir$BNJR) * noise_level),
BNJR = pmax(0, BNJR + Noise)
)
}
# --- SKENARIO 1: Business as Usual (BAU) ---
scen_bau <- data_future_dynamic %>%
mutate(Skenario = "Business as Usual (BAU)") %>%
calc_flood()
# --- SKENARIO 2: Kebijakan Biasa (Moderasi DEF) ---
# Logika: Jika DEF > 0 (rusak), kurangi 50%. Jika DEF <= 0 (reforestasi), biarkan.
scen_modat <- data_future_dynamic %>%
mutate(
DEF = if_else(DEF > 0, DEF * 0.5, DEF),
Skenario = "Kebijakan Biasa (Reduksi DEF 50%)"
) %>%
calc_flood()
# --- SKENARIO 3: Kebijakan Total (Stop DEF & Perbaikan Sungai) ---
# Logika: Jika DEF > 0 (rusak), paksa jadi 0 (Stop). Jika DEF <= 0 (reforestasi), biarkan.
scen_total <- data_future_dynamic %>%
mutate(
DEF = if_else(DEF > 0, 0, DEF),
SNGI = SNGI * 0.8,
Skenario = "Kebijakan Total (Stop DEF + Drainase)"
) %>%
calc_flood()
# --- Menggabungkan Semua Skenario ---
data_scen_all <- rbind(scen_bau, scen_modat, scen_total)
#===============================================================================
# 13. PERSIAPAN: TITIK SAMBUNG (ANCHOR) AGAR GARIS TIDAK TERPUTUS
#===============================================================================
# Mengambil data tahun 2024 sebagai titik awal untuk semua proyeksi
data_anchor_2024 <- banjir %>%
select(PROV, TAHUN, BNJR) %>%
filter(TAHUN == 2024) %>%
mutate(TAHUN = as.numeric(as.character(TAHUN)))
# Menyiapkan data tiap skenario agar menyambung dari 2024
data_bau_conn <- rbind(data_anchor_2024, scen_bau %>% select(PROV, TAHUN, BNJR)) %>%
mutate(Skenario = "Business as Usual (BAU)")
data_modat_conn <- rbind(data_anchor_2024, scen_modat %>% select(PROV, TAHUN, BNJR)) %>%
mutate(Skenario = "Kebijakan Biasa (Reduksi DEF 50%)")
data_total_conn <- rbind(data_anchor_2024, scen_total %>% select(PROV, TAHUN, BNJR)) %>%
mutate(Skenario = "Kebijakan Total (Stop DEF + Drainase)")
data_hist_plot <- banjir %>%
select(PROV, TAHUN, BNJR) %>%
mutate(TAHUN = as.numeric(as.character(TAHUN)), Skenario = "Historis")
#===============================================================================
# 14. VISUALISASI PERBANDINGAN SKENARIO (RATA-RATA SUMATERA)
#===============================================================================
# Menghitung rata-rata per skenario untuk tren umum
avg_hist <- data_hist_plot %>% group_by(TAHUN) %>% summarise(BNJR = mean(BNJR)) %>% mutate(Skenario = "Historis")
avg_bau <- data_bau_conn %>% group_by(TAHUN) %>% summarise(BNJR = mean(BNJR)) %>% mutate(Skenario = "Business as Usual (BAU)")
avg_modat <- data_modat_conn %>% group_by(TAHUN) %>% summarise(BNJR = mean(BNJR)) %>% mutate(Skenario = "Kebijakan Biasa (Reduksi DEF 50%)")
avg_total <- data_total_conn %>% group_by(TAHUN) %>% summarise(BNJR = mean(BNJR)) %>% mutate(Skenario = "Kebijakan Total (Stop DEF + Drainase)")
scen_avg_connected <- rbind(avg_hist, avg_bau, avg_modat, avg_total)
p_scen <- ggplot(scen_avg_connected, aes(x = TAHUN, y = BNJR, color = Skenario, group = Skenario)) +
geom_line(aes(linetype = Skenario), linewidth = 1.2) +
geom_point(size = 2) +
scale_x_continuous(breaks = 2018:2030) +
scale_color_manual(values = c("Historis" = "black", "Business as Usual (BAU)" = "#e41a1c",
"Kebijakan Biasa (Reduksi DEF 50%)" = "#377eb8",
"Kebijakan Total (Stop DEF + Drainase)" = "#4daf4a")) +
scale_linetype_manual(values = c("Historis" = "solid", "Business as Usual (BAU)" = "dashed",
"Kebijakan Biasa (Reduksi DEF 50%)" = "dashed",
"Kebijakan Total (Stop DEF + Drainase)" = "dashed")) +
labs(title = "Rata-rata Dampak Kebijakan terhadap Banjir di Sumatera",
subtitle = "Garis menyambung dari data historis ke simulasi 2025-2030",
x = "Tahun", y = "Rata-rata Kejadian Banjir") +
theme_minimal() + theme(legend.position = "bottom")
print(p_scen)
#===============================================================================
# 15. VISUALISASI PERBANDINGAN SKENARIO (FORMAT TIME SERIES PER PROVINSI)
#===============================================================================
# Menggabungkan semua data per provinsi
data_scen_viz_detail <- rbind(data_hist_plot, data_bau_conn, data_modat_conn, data_total_conn)
p_scen_detail <- ggplot(data_scen_viz_detail, aes(x = TAHUN, y = BNJR, color = Skenario, group = Skenario)) +
annotate("rect", xmin = 2024, xmax = 2030, ymin = -Inf, ymax = Inf, fill = "gray95", alpha = 0.5) +
geom_line(aes(linetype = Skenario), linewidth = 0.8) +
geom_point(size = 1, alpha = 0.4) +
facet_wrap(~PROV, scales = "free_y", ncol = 3) +
scale_x_continuous(breaks = c(2018, 2024, 2030)) +
scale_color_manual(values = c("Historis" = "black", "Business as Usual (BAU)" = "#e41a1c",
"Kebijakan Biasa (Reduksi DEF 50%)" = "#377eb8",
"Kebijakan Total (Stop DEF + Drainase)" = "#4daf4a")) +
scale_linetype_manual(values = c("Historis" = "solid", "Business as Usual (BAU)" = "dashed",
"Kebijakan Biasa (Reduksi DEF 50%)" = "dashed",
"Kebijakan Total (Stop DEF + Drainase)" = "dashed")) +
labs(title = "Detail Skenario Banjir per Provinsi di Sumatera (2018-2030)",
subtitle = "Analisis per Wilayah: Efektivitas Intervensi Kebijakan",
x = "Tahun", y = "Jumlah Kejadian Banjir") +
theme_minimal() +
theme(legend.position = "bottom", strip.background = element_rect(fill = "gray20"),
strip.text = element_text(color = "white", face = "bold"),
axis.text.x = element_text(size = 8))
print(p_scen_detail)
#===============================================================================
# 16. VISUALISASI KUADRAN SPASIAL (2018, 2024, 2030 BAU, 2030 TOTAL)
#===============================================================================
# 1. Menyiapkan 4 Dataset untuk 4 Kuadran
# A. Akhir Historis (2024)
data_2024 <- gwpr_summary_sig %>% filter(TAHUN == 2024) %>%
select(PROV, BNJR = y, DEF, LAT, LONG) %>% mutate(LABEL = "2024: Kondisi Eksisting")
# B. Proyeksi 2030 Tanpa Kebijakan (BAU)
data_2030_bau <- scen_bau %>% filter(TAHUN == 2030) %>%
left_join(wilayah_raw %>% select(name, latitude, longitude), by = c("PROV" = "name")) %>%
select(PROV, BNJR, DEF, LAT = latitude, LONG = longitude) %>%
mutate(LABEL = "2030: Tanpa Kebijakan (BAU)")
# C. Proyeksi 2030 Tanpa Kebijakan (BAU)
data_2030_biasa <- scen_modat %>% filter(TAHUN == 2030) %>%
left_join(wilayah_raw %>% select(name, latitude, longitude), by = c("PROV" = "name")) %>%
select(PROV, BNJR, DEF, LAT = latitude, LONG = longitude) %>%
mutate(LABEL = "2030: Kebijakan")
# D. Proyeksi 2030 Dengan Kebijakan Total
data_2030_total <- scen_total %>% filter(TAHUN == 2030) %>%
left_join(wilayah_raw %>% select(name, latitude, longitude), by = c("PROV" = "name")) %>%
select(PROV, BNJR, DEF, LAT = latitude, LONG = longitude) %>%
mutate(LABEL = "2030: Kebijakan Total (Stop DEF)")
# 2. Gabungkan & Konversi ke SF
peta_kuadran_data <- rbind(data_2024, data_2030_bau, data_2030_biasa, data_2030_total)
# Mengatur urutan agar 2018 & 2024 di atas, 2030 di bawah
peta_kuadran_data$LABEL <- factor(peta_kuadran_data$LABEL,
levels = c("2024: Kondisi Eksisting", "2030: Tanpa Kebijakan (BAU)",
"2030: Kebijakan", "2030: Kebijakan Total (Stop DEF)"))
peta_kuadran_sf <- st_as_sf(peta_kuadran_data, coords = c("LONG", "LAT"), crs = 4326)
# 3. Plotting 4 Kuadran
p_perbandingan_final <- ggplot() +
geom_sf(data = sumatera_map, fill = "gray97", color = "gray85", size = 0.1) +
geom_sf(data = peta_kuadran_sf, aes(color = BNJR, size = DEF), alpha = 0.8) +
# Facet wrap 2x2
facet_wrap(~LABEL, ncol = 2) +
scale_color_viridis_c(option = "plasma", name = "Intensitas Banjir") +
scale_size_continuous(range = c(4, 12), name = "Tingkat Deforestasi") +
labs(title = "Evolusi & Intervensi Spasial Banjir Sumatera (2018-2030)",
subtitle = "Perbandingan Realita Historis vs Dampak Implementasi Kebijakan Total",
caption = "Analisis GWPR: Ukuran titik (Deforestasi) | Warna titik (Frekuensi Banjir)") +
theme_minimal() +
theme(
axis.text = element_blank(),
panel.grid = element_blank(),
strip.background = element_rect(fill = "gray20"),
strip.text = element_text(color = "white", face = "bold", size = 10),
legend.position = "right"
)
print(p_perbandingan_final)
#===============================================================================
# 17. TABEL RINGKASAN HASIL PREDIKSI & ANALISIS DAMPAK KEBIJAKAN (2030)
#===============================================================================
# --- 17.1. Tabel Efektivitas Skenario Secara Agregat (Sumatera) ---
tabel_efektivitas <- data_scen_all %>%
filter(TAHUN == 2030) %>%
group_by(Skenario) %>%
summarise(
Rerata_Banjir = mean(BNJR),
Min_Banjir = min(BNJR),
Max_Banjir = max(BNJR)
) %>%
mutate(
Persen_Penurunan = (1 - Rerata_Banjir / Rerata_Banjir[Skenario == "Business as Usual (BAU)"]) * 100,
Status = c("Kritis (Baseline)", "Waspada", "Terkendali") # Urutan sesuai nama skenario
)
kable(tabel_efektivitas, digits = 2, format = "markdown",
caption = "Ringkasan Prediksi Banjir Sumatera Tahun 2030 Berdasarkan Skenario")
Ringkasan Prediksi Banjir Sumatera Tahun 2030 Berdasarkan
Skenario
| Business as Usual (BAU) |
54.22 |
0 |
232.00 |
0.00 |
Kritis (Baseline) |
| Kebijakan Biasa (Reduksi DEF 50%) |
46.03 |
0 |
220.21 |
15.11 |
Waspada |
| Kebijakan Total (Stop DEF + Drainase) |
37.34 |
0 |
204.56 |
31.13 |
Terkendali |
# --- 17.2. Tabel Dampak Spesifik per Provinsi (Business as Usual vs Kebijakan Total) ---
# Membandingkan kondisi 2024 (Eksisting) dengan 2030 (BAU dan Total)
perbandingan_provinsi <- data_scen_all %>%
filter(TAHUN == 2030) %>%
select(PROV, Skenario, BNJR) %>%
tidyr::pivot_wider(names_from = Skenario, values_from = BNJR) %>%
mutate(
Selisih_Penyelamatan = `Business as Usual (BAU)` - `Kebijakan Total (Stop DEF + Drainase)`,
Persen_Efektivitas = (Selisih_Penyelamatan / `Business as Usual (BAU)`) * 100
) %>%
arrange(desc(Selisih_Penyelamatan))
kable(perbandingan_provinsi, digits = 2, format = "markdown",
caption = "Analisis Efektivitas Kebijakan Total terhadap Penyelamatan Wilayah (Tahun 2030)")
Analisis Efektivitas Kebijakan Total terhadap Penyelamatan
Wilayah (Tahun 2030)
| ACEH |
43.44 |
14.45 |
0.00 |
43.44 |
100.00 |
| SUMATERA BARAT |
108.15 |
91.08 |
71.48 |
36.67 |
33.91 |
| SUMATERA UTARA |
232.00 |
220.21 |
204.56 |
27.44 |
11.83 |
| KEPULAUAN BANGKA BELITUNG |
21.82 |
4.16 |
0.00 |
21.82 |
100.00 |
| BENGKULU |
36.93 |
32.73 |
22.48 |
14.44 |
39.12 |
| LAMPUNG |
84.02 |
84.02 |
71.49 |
12.53 |
14.91 |
| JAMBI |
7.72 |
7.72 |
0.51 |
7.21 |
93.44 |
| KEPULAUAN RIAU |
8.16 |
5.92 |
2.92 |
5.24 |
64.16 |
| RIAU |
0.00 |
0.00 |
0.00 |
0.00 |
NaN |
| SUMATERA SELATAN |
0.00 |
0.00 |
0.00 |
0.00 |
NaN |
# --- 17.3. Tabel Tren Sensitivitas Koefisien Lokal (Output GWPR Dinamis) ---
# Mengambil rata-rata koefisien prediksi untuk melihat variabel mana yang semakin berbahaya
tabel_tren_koef <- data_future_dynamic %>%
filter(TAHUN %in% c(2025, 2030)) %>%
group_by(TAHUN) %>%
summarise(
Pengaruh_DEF = mean(DEF_pred),
Pengaruh_KPDT = mean(KPDT_pred),
Pengaruh_SNGI = mean(SNGI_pred)
)
kable(tabel_tren_koef, digits = 4, format = "markdown",
caption = "Evolusi Sensitivitas (Koefisien) Pengaruh Variabel terhadap Banjir")
Evolusi Sensitivitas (Koefisien) Pengaruh Variabel terhadap
Banjir
| 2025 |
0.0037 |
-0.1602 |
0.0117 |
| 2030 |
0.0069 |
-0.2693 |
0.0189 |
# --- 17.4. Tabel Akumulasi Beban Banjir (Total Kejadian 2025-2030) ---
# Menghitung total perkiraan banjir yang akan terjadi selama 6 tahun ke depan
tabel_kumulatif <- data_scen_all %>%
group_by(Skenario) %>%
summarise(
Total_Kejadian_6Thn = sum(BNJR),
Rata_Rata_Tahunan = mean(BNJR),
Beban_Maksimal_Prov = max(BNJR)
) %>%
mutate(Potensi_Nyawa_Terselamatkan = paste0(round((1 - Total_Kejadian_6Thn / Total_Kejadian_6Thn[1])*100, 1), "%"))
kable(tabel_kumulatif, digits = 2, format = "markdown",
caption = "Proyeksi Akumulasi Beban Kejadian Banjir di Sumatera (2025-2030)")
Proyeksi Akumulasi Beban Kejadian Banjir di Sumatera
(2025-2030)
| Business as Usual (BAU) |
3244.06 |
54.07 |
232.00 |
0% |
| Kebijakan Biasa (Reduksi DEF 50%) |
2846.79 |
47.45 |
220.21 |
12.2% |
| Kebijakan Total (Stop DEF + Drainase) |
2316.26 |
38.60 |
204.56 |
28.6% |
# --- 17.5. Identifikasi 'Top 5 Hotspot' Prioritas (Skenario BAU 2030) ---
# Provinsi mana yang paling kritis jika tidak ada intervensi sama sekali
hotspot_2030 <- scen_bau %>%
filter(TAHUN == 2030) %>%
arrange(desc(BNJR)) %>%
slice(1:5) %>%
select(PROV, Prediksi_Banjir_2030 = BNJR, Deforestasi_Proyeksi = DEF)
kable(hotspot_2030, digits = 2, format = "markdown",
caption = "Top 5 Provinsi dengan Risiko Banjir Tertinggi pada Tahun 2030 (Skenario BAU)")
Top 5 Provinsi dengan Risiko Banjir Tertinggi pada Tahun 2030
(Skenario BAU)
| SUMATERA UTARA |
232.00 |
3456.01 |
| SUMATERA BARAT |
108.15 |
5233.22 |
| LAMPUNG |
84.02 |
-931.75 |
| ACEH |
43.44 |
8348.52 |
| BENGKULU |
36.93 |
1190.14 |
# --- 17.6. Matriks Kontribusi Variabel terhadap Kenaikan Banjir ---
# Mengukur variabel mana yang paling 'jahat' (sensitif) di masa depan
matriks_sensitivitas <- data_future_dynamic %>%
summarise(
Kontribusi_DEF = sum(abs(DEF * DEF_pred)),
Kontribusi_KPDT = sum(abs(KPDT * KPDT_pred)),
Kontribusi_Drainase = sum(abs(SNGI * SNGI_pred))
) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Faktor_Pemicu", values_to = "Nilai_Kontribusi") %>%
mutate(Persentase_Andil = (Nilai_Kontribusi / sum(Nilai_Kontribusi)) * 100) %>%
arrange(desc(Persentase_Andil))
kable(matriks_sensitivitas, digits = 2, format = "markdown",
caption = "Matriks Faktor Pemicu Utama Banjir di Sumatera (Proyeksi 2025-2030)")
Matriks Faktor Pemicu Utama Banjir di Sumatera (Proyeksi
2025-2030)
| Kontribusi_DEF |
2666.58 |
44.21 |
| Kontribusi_KPDT |
1942.51 |
32.21 |
| Kontribusi_Drainase |
1421.93 |
23.58 |
# --- 17.4. Tabel Akumulasi Beban Banjir (Total Kejadian 2025-2030) ---
# Menghitung total perkiraan banjir yang akan terjadi selama 6 tahun ke depan
tabel_kumulatif <- data_scen_all %>%
group_by(Skenario) %>%
summarise(
Total_Kejadian_6Thn = sum(BNJR),
Rata_Rata_Tahunan = mean(BNJR),
Beban_Maksimal_Prov = max(BNJR)
) %>%
mutate(Potensi_Nyawa_Terselamatkan = paste0(round((1 - Total_Kejadian_6Thn / Total_Kejadian_6Thn[1])*100, 1), "%"))
kable(tabel_kumulatif, digits = 2, format = "markdown",
caption = "Proyeksi Akumulasi Beban Kejadian Banjir di Sumatera (2025-2030)")
Proyeksi Akumulasi Beban Kejadian Banjir di Sumatera
(2025-2030)
| Business as Usual (BAU) |
3244.06 |
54.07 |
232.00 |
0% |
| Kebijakan Biasa (Reduksi DEF 50%) |
2846.79 |
47.45 |
220.21 |
12.2% |
| Kebijakan Total (Stop DEF + Drainase) |
2316.26 |
38.60 |
204.56 |
28.6% |
##===============================================================================
# 17.5. IDENTIFIKASI HOTSPOT PRIORITAS & PERBANDINGAN ANTAR SKENARIO (2030)
#===============================================================================
##===============================================================================
# 17.5. TABEL LENGKAP: PERBANDINGAN BANJIR & DEFORESTASI (3 SKENARIO 2030)
#===============================================================================
# 1. Ambil Top 5 Provinsi Risiko Tertinggi
top_5_prov <- scen_bau %>%
filter(TAHUN == 2030) %>%
arrange(desc(BNJR)) %>%
slice(1:5) %>%
pull(PROV)
# 2. Gabungkan data dengan semua kolom Deforestasi (BAU, Biasa, Total)
hotspot_lengkap <- data_scen_all %>%
filter(TAHUN == 2030, PROV %in% top_5_prov) %>%
select(PROV, Skenario, BNJR, DEF) %>%
tidyr::pivot_wider(
names_from = Skenario,
values_from = c(BNJR, DEF),
names_glue = "{Skenario}_{.value}"
) %>%
# Mengatur ulang dan memberi nama kolom yang informatif
select(
Provinsi = PROV,
`Bjr_BAU` = `Business as Usual (BAU)_BNJR`,
`Bjr_Biasa` = `Kebijakan Biasa (Reduksi DEF 50%)_BNJR`,
`Bjr_Total` = `Kebijakan Total (Stop DEF + Drainase)_BNJR`,
`Def_BAU` = `Business as Usual (BAU)_DEF`,
`Def_Biasa` = `Kebijakan Biasa (Reduksi DEF 50%)_DEF`, # Ini yang Anda cari
`Def_Total` = `Kebijakan Total (Stop DEF + Drainase)_DEF`
) %>%
arrange(desc(Bjr_BAU))
# 3. Tampilkan dengan kable
kable(hotspot_lengkap, digits = 2, format = "markdown",
caption = "Tabel Perbandingan Dampak Deforestasi terhadap Kejadian Banjir Tahun 2030")
Tabel Perbandingan Dampak Deforestasi terhadap Kejadian Banjir
Tahun 2030
| SUMATERA UTARA |
232.00 |
220.21 |
204.56 |
3456.01 |
1728.01 |
0.00 |
| SUMATERA BARAT |
108.15 |
91.08 |
71.48 |
5233.22 |
2616.61 |
0.00 |
| LAMPUNG |
84.02 |
84.02 |
71.49 |
-931.75 |
-931.75 |
-931.75 |
| ACEH |
43.44 |
14.45 |
0.00 |
8348.52 |
4174.26 |
0.00 |
| BENGKULU |
36.93 |
32.73 |
22.48 |
1190.14 |
595.07 |
0.00 |
# --- 17.6. Matriks Kontribusi Variabel terhadap Kenaikan Banjir ---
# Mengukur variabel mana yang paling 'jahat' (sensitif) di masa depan
matriks_sensitivitas <- data_future_dynamic %>%
summarise(
Kontribusi_DEF = sum(abs(DEF * DEF_pred)),
Kontribusi_KPDT = sum(abs(KPDT * KPDT_pred)),
Kontribusi_Drainase = sum(abs(SNGI * SNGI_pred))
) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Faktor_Pemicu", values_to = "Nilai_Kontribusi") %>%
mutate(Persentase_Andil = (Nilai_Kontribusi / sum(Nilai_Kontribusi)) * 100) %>%
arrange(desc(Persentase_Andil))
kable(matriks_sensitivitas, digits = 2, format = "markdown",
caption = "Matriks Faktor Pemicu Utama Banjir di Sumatera (Proyeksi 2025-2030)")
Matriks Faktor Pemicu Utama Banjir di Sumatera (Proyeksi
2025-2030)
| Kontribusi_DEF |
2666.58 |
44.21 |
| Kontribusi_KPDT |
1942.51 |
32.21 |
| Kontribusi_Drainase |
1421.93 |
23.58 |