#===============================================================================
# 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_final.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.075507 1.004071 1.277543 1.354044
# 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 |
1.3783 |
0.2404 |
Tidak Signifikan |
| Time |
4.2342 |
0.0396 |
Signifikan |
| Twoways |
5.6125 |
0.0604 |
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 |
14979.09 |
604.2655 |
637.9929 |
274.1363 |
0.3180 |
| Time |
17704.37 |
609.9664 |
636.9484 |
285.5465 |
0.6991 |
| Twoways |
13817.47 |
610.6150 |
657.8334 |
282.3034 |
0.1089 |
# 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.53333 |
0.02998 |
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 |
1.3783 |
0.2404 |
Tidak Signifikan |
| Time |
4.2342 |
0.0396 |
Signifikan |
| Twoways |
5.6125 |
0.0604 |
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.0665 |
0.9995 |
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.0973 |
0.4909 |
Terpenuhi |
| Autokorelasi (Breusch-Godfrey) |
11.0656 |
0.1358 |
Terpenuhi |
| Heteroskedastisitas (Breusch-Pagan) |
1.3783 |
0.2404 |
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: 29786.28
## Adaptive bandwidth: 39 CV score: 24709.97
## Adaptive bandwidth: 30 CV score: 24815.81
## Adaptive bandwidth: 42 CV score: 24709.97
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: 27114.28
## Adaptive bandwidth: 39 CV score: 27675.33
## Adaptive bandwidth: 58 CV score: 27045.36
## Adaptive bandwidth: 62 CV score: 27045.36
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: 27393.88
## Adaptive bandwidth: 39 CV score: 27558.25
## Adaptive bandwidth: 58 CV score: 27316.14
## Adaptive bandwidth: 62 CV score: 27316.14
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: 29905.5
## Fixed bandwidth: 4.867243 CV score: 27727.02
## Fixed bandwidth: 3.009095 CV score: 109660.3
## Fixed bandwidth: 6.015642 CV score: 31593.56
## Fixed bandwidth: 4.157494 CV score: 28036.67
## Fixed bandwidth: 5.305893 CV score: 30253.19
## Fixed bandwidth: 4.596143 CV score: 27081.19
## Fixed bandwidth: 4.428594 CV score: 27242.9
## Fixed bandwidth: 4.699694 CV score: 27083.14
## Fixed bandwidth: 4.532145 CV score: 27121.17
## Fixed bandwidth: 4.635696 CV score: 27067.98
## Fixed bandwidth: 4.660141 CV score: 27063.8
## Fixed bandwidth: 4.675249 CV score: 27065.83
## Fixed bandwidth: 4.650804 CV score: 27065.03
## Fixed bandwidth: 4.665912 CV score: 27063.84
## Fixed bandwidth: 4.656575 CV score: 27064.19
## Fixed bandwidth: 4.662345 CV score: 27063.71
## Fixed bandwidth: 4.663708 CV score: 27063.72
## Fixed bandwidth: 4.661503 CV score: 27063.73
## Fixed bandwidth: 4.662866 CV score: 27063.7
## Fixed bandwidth: 4.663187 CV score: 27063.71
## Fixed bandwidth: 4.662667 CV score: 27063.7
## Fixed bandwidth: 4.662989 CV score: 27063.7
## Fixed bandwidth: 4.66279 CV score: 27063.7
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: 26943.29
## Fixed bandwidth: 4.867243 CV score: 27682.46
## Fixed bandwidth: 9.731939 CV score: 26796.68
## Fixed bandwidth: 10.88034 CV score: 26743.31
## Fixed bandwidth: 11.59009 CV score: 26718.46
## Fixed bandwidth: 12.02874 CV score: 26705.39
## Fixed bandwidth: 12.29984 CV score: 26698.05
## Fixed bandwidth: 12.46739 CV score: 26693.77
## Fixed bandwidth: 12.57094 CV score: 26691.22
## Fixed bandwidth: 12.63494 CV score: 26689.67
## Fixed bandwidth: 12.67449 CV score: 26688.73
## Fixed bandwidth: 12.69893 CV score: 26688.15
## Fixed bandwidth: 12.71404 CV score: 26687.79
## Fixed bandwidth: 12.72338 CV score: 26687.57
## Fixed bandwidth: 12.72915 CV score: 26687.44
## Fixed bandwidth: 12.73272 CV score: 26687.35
## Fixed bandwidth: 12.73492 CV score: 26687.3
## Fixed bandwidth: 12.73628 CV score: 26687.27
## Fixed bandwidth: 12.73712 CV score: 26687.25
## Fixed bandwidth: 12.73764 CV score: 26687.24
## Fixed bandwidth: 12.73797 CV score: 26687.23
## Fixed bandwidth: 12.73816 CV score: 26687.23
## Fixed bandwidth: 12.73829 CV score: 26687.22
## Fixed bandwidth: 12.73836 CV score: 26687.22
## Fixed bandwidth: 12.73841 CV score: 26687.22
## Fixed bandwidth: 12.73844 CV score: 26687.22
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: 27318.37
## Fixed bandwidth: 4.867243 CV score: 27700.34
## Fixed bandwidth: 9.731939 CV score: 27175.69
## Fixed bandwidth: 10.88034 CV score: 27109.75
## Fixed bandwidth: 11.59009 CV score: 27075.16
## Fixed bandwidth: 12.02874 CV score: 27055.74
## Fixed bandwidth: 12.29984 CV score: 27044.4
## Fixed bandwidth: 12.46739 CV score: 27037.63
## Fixed bandwidth: 12.57094 CV score: 27033.54
## Fixed bandwidth: 12.63494 CV score: 27031.04
## Fixed bandwidth: 12.67449 CV score: 27029.51
## Fixed bandwidth: 12.69893 CV score: 27028.57
## Fixed bandwidth: 12.71404 CV score: 27027.99
## Fixed bandwidth: 12.72338 CV score: 27027.63
## Fixed bandwidth: 12.72915 CV score: 27027.41
## Fixed bandwidth: 12.73272 CV score: 27027.27
## Fixed bandwidth: 12.73492 CV score: 27027.19
## Fixed bandwidth: 12.73628 CV score: 27027.13
## Fixed bandwidth: 12.73712 CV score: 27027.1
## Fixed bandwidth: 12.73764 CV score: 27027.08
## Fixed bandwidth: 12.73797 CV score: 27027.07
## Fixed bandwidth: 12.73816 CV score: 27027.06
## Fixed bandwidth: 12.73829 CV score: 27027.06
## Fixed bandwidth: 12.73836 CV score: 27027.05
## Fixed bandwidth: 12.73841 CV score: 27027.05
## Fixed bandwidth: 12.73844 CV score: 27027.05
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 |
16908.11 |
0.741 |
0.658 |
597.621 |
624.584 |
- |
| Adaptive Gaussian |
21915.45 |
0.664 |
0.631 |
606.601 |
615.954 |
- |
| Adaptive Exponential |
21284.78 |
0.674 |
0.631 |
605.631 |
616.631 |
- |
| Fixed Bisquare |
16491.87 |
0.747 |
0.668 |
595.886 |
622.872 |
Model Terbaik |
| Fixed Gaussian |
22129.70 |
0.661 |
0.632 |
606.858 |
615.588 |
- |
| Fixed Exponential |
21725.95 |
0.667 |
0.630 |
606.382 |
616.320 |
- |
#===============================================================================
# 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 = "bisquare",
adaptive = FALSE)
# 4. Estimasi Model GWR
model_gwr <- gwr.basic(BNJR ~ DEF + KPDT + CAC + SNGI,
data = data_tahun_ini,
bw = bw_optimum,
kernel = "bisquare",
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: 3448.093
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 2445.578
## Fixed bandwidth: 10.88034 CV score: 2342.993
## Fixed bandwidth: 11.59009 CV score: 2494.127
## Fixed bandwidth: 10.44169 CV score: 2283.551
## Fixed bandwidth: 10.17059 CV score: 2296.534
## Fixed bandwidth: 10.60924 CV score: 2298.521
## Fixed bandwidth: 10.33814 CV score: 2281.973
## Fixed bandwidth: 10.27414 CV score: 2284.765
## Fixed bandwidth: 10.37769 CV score: 2281.75
## Fixed bandwidth: 10.40214 CV score: 2282.138
## Fixed bandwidth: 10.36258 CV score: 2281.708
## Fixed bandwidth: 10.35325 CV score: 2281.76
## Fixed bandwidth: 10.36835 CV score: 2281.706
## Fixed bandwidth: 10.37192 CV score: 2281.716
## Fixed bandwidth: 10.36615 CV score: 2281.704
## Fixed bandwidth: 10.36479 CV score: 2281.705
## Fixed bandwidth: 10.36699 CV score: 2281.705
## Fixed bandwidth: 10.36563 CV score: 2281.704
## Selesai memproses tahun: 2018
## Fixed bandwidth: 7.873791 CV score: 4738.214
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 2541.909
## Fixed bandwidth: 10.88034 CV score: 2313.338
## Fixed bandwidth: 11.59009 CV score: 2203.362
## Fixed bandwidth: 12.02874 CV score: 2149.135
## Fixed bandwidth: 12.29984 CV score: 2120.473
## Fixed bandwidth: 12.46739 CV score: 2105.409
## Fixed bandwidth: 12.57094 CV score: 2097.042
## Fixed bandwidth: 12.63494 CV score: 2092.176
## Fixed bandwidth: 12.67449 CV score: 2089.275
## Fixed bandwidth: 12.69893 CV score: 2087.519
## Fixed bandwidth: 12.71404 CV score: 2086.448
## Fixed bandwidth: 12.72338 CV score: 2085.791
## Fixed bandwidth: 12.72915 CV score: 2085.388
## Fixed bandwidth: 12.73272 CV score: 2085.139
## Fixed bandwidth: 12.73492 CV score: 2084.985
## Fixed bandwidth: 12.73628 CV score: 2084.89
## Fixed bandwidth: 12.73712 CV score: 2084.832
## Fixed bandwidth: 12.73764 CV score: 2084.796
## Fixed bandwidth: 12.73797 CV score: 2084.773
## Fixed bandwidth: 12.73816 CV score: 2084.76
## Fixed bandwidth: 12.73829 CV score: 2084.751
## Fixed bandwidth: 12.73836 CV score: 2084.746
## Fixed bandwidth: 12.73841 CV score: 2084.743
## Fixed bandwidth: 12.73844 CV score: 2084.741
## Selesai memproses tahun: 2019
## Fixed bandwidth: 7.873791 CV score: 240428.7
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 12833.07
## Fixed bandwidth: 10.88034 CV score: 5977.298
## Fixed bandwidth: 11.59009 CV score: 5230.679
## Fixed bandwidth: 12.02874 CV score: 5005.014
## Fixed bandwidth: 12.29984 CV score: 4906.422
## Fixed bandwidth: 12.46739 CV score: 4849.007
## Fixed bandwidth: 12.57094 CV score: 4814.231
## Fixed bandwidth: 12.63494 CV score: 4793.183
## Fixed bandwidth: 12.67449 CV score: 4780.378
## Fixed bandwidth: 12.69893 CV score: 4772.547
## Fixed bandwidth: 12.71404 CV score: 4767.741
## Fixed bandwidth: 12.72338 CV score: 4764.784
## Fixed bandwidth: 12.72915 CV score: 4762.961
## Fixed bandwidth: 12.73272 CV score: 4761.837
## Fixed bandwidth: 12.73492 CV score: 4761.143
## Fixed bandwidth: 12.73628 CV score: 4760.714
## Fixed bandwidth: 12.73712 CV score: 4760.449
## Fixed bandwidth: 12.73764 CV score: 4760.285
## Fixed bandwidth: 12.73797 CV score: 4760.184
## Fixed bandwidth: 12.73816 CV score: 4760.122
## Fixed bandwidth: 12.73829 CV score: 4760.083
## Fixed bandwidth: 12.73836 CV score: 4760.059
## Fixed bandwidth: 12.73841 CV score: 4760.045
## Fixed bandwidth: 12.73844 CV score: 4760.035
## Selesai memproses tahun: 2020
## Fixed bandwidth: 7.873791 CV score: 13100.85
## Fixed bandwidth: 4.867243 CV score: 177247.7
## Fixed bandwidth: 9.731939 CV score: 2808.59
## Fixed bandwidth: 10.88034 CV score: 2409.504
## Fixed bandwidth: 11.59009 CV score: 2251.613
## Fixed bandwidth: 12.02874 CV score: 2180.191
## Fixed bandwidth: 12.29984 CV score: 2143.658
## Fixed bandwidth: 12.46739 CV score: 2123.789
## Fixed bandwidth: 12.57094 CV score: 2112.462
## Fixed bandwidth: 12.63494 CV score: 2105.796
## Fixed bandwidth: 12.67449 CV score: 2101.799
## Fixed bandwidth: 12.69893 CV score: 2099.373
## Fixed bandwidth: 12.71404 CV score: 2097.891
## Fixed bandwidth: 12.72338 CV score: 2096.981
## Fixed bandwidth: 12.72915 CV score: 2096.421
## Fixed bandwidth: 12.73272 CV score: 2096.076
## Fixed bandwidth: 12.73492 CV score: 2095.863
## Fixed bandwidth: 12.73628 CV score: 2095.732
## Fixed bandwidth: 12.73712 CV score: 2095.651
## Fixed bandwidth: 12.73764 CV score: 2095.6
## Fixed bandwidth: 12.73797 CV score: 2095.569
## Fixed bandwidth: 12.73816 CV score: 2095.55
## Fixed bandwidth: 12.73829 CV score: 2095.538
## Fixed bandwidth: 12.73836 CV score: 2095.531
## Fixed bandwidth: 12.73841 CV score: 2095.527
## Fixed bandwidth: 12.73844 CV score: 2095.524
## Selesai memproses tahun: 2021
## Fixed bandwidth: 7.873791 CV score: 34317.57
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 10152.61
## Fixed bandwidth: 10.88034 CV score: 7846.501
## Fixed bandwidth: 11.59009 CV score: 7165.265
## Fixed bandwidth: 12.02874 CV score: 6902.448
## Fixed bandwidth: 12.29984 CV score: 6779.457
## Fixed bandwidth: 12.46739 CV score: 6716.626
## Fixed bandwidth: 12.57094 CV score: 6682.181
## Fixed bandwidth: 12.63494 CV score: 6662.374
## Fixed bandwidth: 12.67449 CV score: 6650.657
## Fixed bandwidth: 12.69893 CV score: 6643.607
## Fixed bandwidth: 12.71404 CV score: 6639.321
## Fixed bandwidth: 12.72338 CV score: 6636.699
## Fixed bandwidth: 12.72915 CV score: 6635.089
## Fixed bandwidth: 12.73272 CV score: 6634.098
## Fixed bandwidth: 12.73492 CV score: 6633.486
## Fixed bandwidth: 12.73628 CV score: 6633.109
## Fixed bandwidth: 12.73712 CV score: 6632.876
## Fixed bandwidth: 12.73764 CV score: 6632.732
## Fixed bandwidth: 12.73797 CV score: 6632.643
## Fixed bandwidth: 12.73816 CV score: 6632.588
## Fixed bandwidth: 12.73829 CV score: 6632.555
## Fixed bandwidth: 12.73836 CV score: 6632.534
## Fixed bandwidth: 12.73841 CV score: 6632.521
## Fixed bandwidth: 12.73844 CV score: 6632.513
## Selesai memproses tahun: 2022
## Fixed bandwidth: 7.873791 CV score: 54951.4
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 31918.66
## Fixed bandwidth: 10.88034 CV score: 24311.04
## Fixed bandwidth: 11.59009 CV score: 22156.89
## Fixed bandwidth: 12.02874 CV score: 21288.13
## Fixed bandwidth: 12.29984 CV score: 20855.83
## Fixed bandwidth: 12.46739 CV score: 20579.21
## Fixed bandwidth: 12.57094 CV score: 20402.05
## Fixed bandwidth: 12.63494 CV score: 20291.38
## Fixed bandwidth: 12.67449 CV score: 20222.77
## Fixed bandwidth: 12.69893 CV score: 20180.34
## Fixed bandwidth: 12.71404 CV score: 20154.12
## Fixed bandwidth: 12.72338 CV score: 20137.91
## Fixed bandwidth: 12.72915 CV score: 20127.9
## Fixed bandwidth: 12.73272 CV score: 20121.71
## Fixed bandwidth: 12.73492 CV score: 20117.89
## Fixed bandwidth: 12.73628 CV score: 20115.53
## Fixed bandwidth: 12.73712 CV score: 20114.07
## Fixed bandwidth: 12.73764 CV score: 20113.17
## Fixed bandwidth: 12.73797 CV score: 20112.61
## Fixed bandwidth: 12.73816 CV score: 20112.26
## Fixed bandwidth: 12.73829 CV score: 20112.05
## Fixed bandwidth: 12.73836 CV score: 20111.92
## Fixed bandwidth: 12.73841 CV score: 20111.84
## Fixed bandwidth: 12.73844 CV score: 20111.79
## Selesai memproses tahun: 2023
## Fixed bandwidth: 7.873791 CV score: 8649.44
## Fixed bandwidth: 4.867243 CV score: Inf
## Fixed bandwidth: 9.731939 CV score: 24763.34
## Fixed bandwidth: 6.725392 CV score: 104506
## Fixed bandwidth: 8.58354 CV score: 7953.768
## Fixed bandwidth: 9.022189 CV score: 27666.72
## Fixed bandwidth: 8.31244 CV score: 9899.306
## Fixed bandwidth: 8.751089 CV score: 6804.46
## Fixed bandwidth: 8.85464 CV score: 27803.93
## Fixed bandwidth: 8.687091 CV score: 17899.32
## Fixed bandwidth: 8.790642 CV score: 8094.184
## Fixed bandwidth: 8.726644 CV score: 48825.76
## Fixed bandwidth: 8.766197 CV score: 6680.04
## Fixed bandwidth: 8.775534 CV score: 11048.94
## Fixed bandwidth: 8.760426 CV score: 8581.459
## Fixed bandwidth: 8.769764 CV score: 19841.4
## Fixed bandwidth: 8.763993 CV score: 13200.93
## Fixed bandwidth: 8.767559 CV score: 7206.858
## Fixed bandwidth: 8.765355 CV score: 6919.519
## Fixed bandwidth: 8.766717 CV score: 29118.54
## Fixed bandwidth: 8.765876 CV score: 32085.01
## Fixed bandwidth: 8.766396 CV score: 58964.16
## Fixed bandwidth: 8.766074 CV score: 12258.22
## Fixed bandwidth: 8.766273 CV score: 6662.834
## Fixed bandwidth: 8.76632 CV score: 32861.42
## Fixed bandwidth: 8.766244 CV score: 9794.049
## 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
|
10
|
10
|
0
|
|
2020
|
0
|
0
|
10
|
9
|
|
2021
|
10
|
10
|
10
|
10
|
|
2022
|
0
|
0
|
2
|
10
|
|
2023
|
9
|
0
|
0
|
1
|
|
2024
|
0
|
0
|
1
|
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 BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, 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 |
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, 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 BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
| 2021 |
Negatif Signifikan |
- |
ACEH, BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG,
KEPULAUAN RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN,
SUMATERA UTARA |
- |
- |
| 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 |
BENGKULU, JAMBI, KEPULAUAN BANGKA BELITUNG, KEPULAUAN
RIAU, LAMPUNG, RIAU, SUMATERA BARAT, SUMATERA SELATAN, SUMATERA
UTARA |
- |
- |
KEPULAUAN RIAU |
| 2023 |
Negatif Signifikan |
- |
- |
- |
- |
| 2024 |
Positif Signifikan |
- |
- |
SUMATERA UTARA |
- |
| 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) ---
# Menggunakan data proyeksi yang sudah kita buat sebelumnya (tren berlanjut)
scen_bau <- data_future_dynamic %>%
mutate(Skenario = "Business as Usual (BAU)") %>%
calc_flood()
# --- SKENARIO 2: Kebijakan Biasa (Moderasi DEF) ---
# Deforestasi dikurangi 50% dari tren biasanya
scen_modat <- data_future_dynamic %>%
mutate(DEF = DEF * 0.5,
Skenario = "Kebijakan Biasa (Reduksi DEF 50%)") %>%
calc_flood()
# --- SKENARIO 3: Kebijakan Total (Stop DEF & Perbaikan Sungai) ---
# Deforestasi = 0, dan variabel Jarak/Normalisasi Sungai (SNGI) membaik 20%
scen_total <- data_future_dynamic %>%
mutate(DEF = 0,
SNGI = SNGI * 0.8, # Asumsi perbaikan manajemen sungai/drainase
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) |
96.57 |
38.79 |
254.79 |
0.00 |
Kritis (Baseline) |
| Kebijakan Biasa (Reduksi DEF 50%) |
81.68 |
0.00 |
240.49 |
15.43 |
Waspada |
| Kebijakan Total (Stop DEF + Drainase) |
61.29 |
0.00 |
218.48 |
36.53 |
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)
| RIAU |
153.05 |
106.99 |
51.91 |
101.15 |
66.09 |
| ACEH |
84.10 |
67.04 |
38.50 |
45.60 |
54.22 |
| KEPULAUAN RIAU |
40.98 |
0.00 |
0.00 |
40.98 |
100.00 |
| SUMATERA UTARA |
254.79 |
240.49 |
218.48 |
36.31 |
14.25 |
| SUMATERA BARAT |
70.92 |
57.94 |
38.80 |
32.12 |
45.29 |
| BENGKULU |
38.79 |
30.81 |
13.34 |
25.45 |
65.60 |
| JAMBI |
49.53 |
43.32 |
25.96 |
23.57 |
47.58 |
| SUMATERA SELATAN |
92.97 |
94.84 |
73.47 |
19.49 |
20.97 |
| LAMPUNG |
132.58 |
133.08 |
118.25 |
14.33 |
10.81 |
| KEPULAUAN BANGKA BELITUNG |
48.03 |
42.23 |
34.21 |
13.82 |
28.77 |
# --- 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.0019 |
-0.2582 |
0.0177 |
| 2030 |
0.0032 |
-0.4461 |
0.0297 |
# --- 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) |
4748.25 |
79.14 |
254.79 |
0% |
| Kebijakan Biasa (Reduksi DEF 50%) |
4049.31 |
67.49 |
240.49 |
14.7% |
| Kebijakan Total (Stop DEF + Drainase) |
3062.17 |
51.04 |
218.48 |
35.5% |
# --- 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 |
254.79 |
9812.38 |
| RIAU |
153.05 |
31716.36 |
| LAMPUNG |
132.58 |
-299.59 |
| SUMATERA SELATAN |
92.97 |
-1152.89 |
| ACEH |
84.10 |
8955.35 |
# --- 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_KPDT |
3558.04 |
48.13 |
| Kontribusi_Drainase |
2301.59 |
31.13 |
| Kontribusi_DEF |
1533.20 |
20.74 |
# --- 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) |
4748.25 |
79.14 |
254.79 |
0% |
| Kebijakan Biasa (Reduksi DEF 50%) |
4049.31 |
67.49 |
240.49 |
14.7% |
| Kebijakan Total (Stop DEF + Drainase) |
3062.17 |
51.04 |
218.48 |
35.5% |
##===============================================================================
# 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 |
254.79 |
240.49 |
218.48 |
9812.38 |
4906.19 |
0 |
| RIAU |
153.05 |
106.99 |
51.91 |
31716.36 |
15858.18 |
0 |
| LAMPUNG |
132.58 |
133.08 |
118.25 |
-299.59 |
-149.80 |
0 |
| SUMATERA SELATAN |
92.97 |
94.84 |
73.47 |
-1152.89 |
-576.44 |
0 |
| ACEH |
84.10 |
67.04 |
38.50 |
8955.35 |
4477.68 |
0 |
# --- 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_KPDT |
3558.04 |
48.13 |
| Kontribusi_Drainase |
2301.59 |
31.13 |
| Kontribusi_DEF |
1533.20 |
20.74 |