Anggota Kelompok
- Ahmad Helmi Yasir (Ketua)
- Najla Khalisa
- Gusti Luqman Nor Hafizh
- Muhammad Ihsan
- Nur Ridho Arif Billah
1. Pustaka (Library)
Tahap pertama memuat seluruh library yang diperlukan untuk analisis data panel, spasial, dan visualisasi.
# GWPR in Sumatra Island 2010-2024 =========================================
# 1. Library ===============================================================
libs <- c("gplots", "kableExtra", "tidyverse", "plm", "tinytex", "tseries",
"lmtest", "broom", "stargazer", "RColorBrewer", "corrplot", "car",
"readxl", "jsonlite", "GWmodel", "sf", "rnaturalearth", "patchwork",
"viridis","plm","knitr","DescTools","tidyr", "ggplot2","spdep")
lapply(libs, library, character.only = TRUE)## [[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"
##
## [[20]]
## [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"
##
## [[21]]
## [1] "knitr" "viridis" "viridisLite" "patchwork"
## [5] "rnaturalearth" "sf" "GWmodel" "Rcpp"
## [9] "sp" "robustbase" "jsonlite" "readxl"
## [13] "car" "carData" "corrplot" "RColorBrewer"
## [17] "stargazer" "broom" "lmtest" "zoo"
## [21] "tseries" "tinytex" "plm" "lubridate"
## [25] "forcats" "stringr" "dplyr" "purrr"
## [29] "readr" "tidyr" "tibble" "ggplot2"
## [33] "tidyverse" "kableExtra" "gplots" "stats"
## [37] "graphics" "grDevices" "utils" "datasets"
## [41] "methods" "base"
##
## [[22]]
## [1] "DescTools" "knitr" "viridis" "viridisLite"
## [5] "patchwork" "rnaturalearth" "sf" "GWmodel"
## [9] "Rcpp" "sp" "robustbase" "jsonlite"
## [13] "readxl" "car" "carData" "corrplot"
## [17] "RColorBrewer" "stargazer" "broom" "lmtest"
## [21] "zoo" "tseries" "tinytex" "plm"
## [25] "lubridate" "forcats" "stringr" "dplyr"
## [29] "purrr" "readr" "tidyr" "tibble"
## [33] "ggplot2" "tidyverse" "kableExtra" "gplots"
## [37] "stats" "graphics" "grDevices" "utils"
## [41] "datasets" "methods" "base"
##
## [[23]]
## [1] "DescTools" "knitr" "viridis" "viridisLite"
## [5] "patchwork" "rnaturalearth" "sf" "GWmodel"
## [9] "Rcpp" "sp" "robustbase" "jsonlite"
## [13] "readxl" "car" "carData" "corrplot"
## [17] "RColorBrewer" "stargazer" "broom" "lmtest"
## [21] "zoo" "tseries" "tinytex" "plm"
## [25] "lubridate" "forcats" "stringr" "dplyr"
## [29] "purrr" "readr" "tidyr" "tibble"
## [33] "ggplot2" "tidyverse" "kableExtra" "gplots"
## [37] "stats" "graphics" "grDevices" "utils"
## [41] "datasets" "methods" "base"
##
## [[24]]
## [1] "DescTools" "knitr" "viridis" "viridisLite"
## [5] "patchwork" "rnaturalearth" "sf" "GWmodel"
## [9] "Rcpp" "sp" "robustbase" "jsonlite"
## [13] "readxl" "car" "carData" "corrplot"
## [17] "RColorBrewer" "stargazer" "broom" "lmtest"
## [21] "zoo" "tseries" "tinytex" "plm"
## [25] "lubridate" "forcats" "stringr" "dplyr"
## [29] "purrr" "readr" "tidyr" "tibble"
## [33] "ggplot2" "tidyverse" "kableExtra" "gplots"
## [37] "stats" "graphics" "grDevices" "utils"
## [41] "datasets" "methods" "base"
##
## [[25]]
## [1] "spdep" "spData" "DescTools" "knitr"
## [5] "viridis" "viridisLite" "patchwork" "rnaturalearth"
## [9] "sf" "GWmodel" "Rcpp" "sp"
## [13] "robustbase" "jsonlite" "readxl" "car"
## [17] "carData" "corrplot" "RColorBrewer" "stargazer"
## [21] "broom" "lmtest" "zoo" "tseries"
## [25] "tinytex" "plm" "lubridate" "forcats"
## [29] "stringr" "dplyr" "purrr" "readr"
## [33] "tidyr" "tibble" "ggplot2" "tidyverse"
## [37] "kableExtra" "gplots" "stats" "graphics"
## [41] "grDevices" "utils" "datasets" "methods"
## [45] "base"2. Persiapan Data
Memuat data kejadian banjir di Pulau Sumatera, melakukan penyesuaian nama kolom, pengaturan tipe data, dan pengecekan missing value, dan mencek outlier dan melakukan penanganan.
# 2. Data ==================================================================
setwd("D:/Kuliah/Lomba/PKM-AI 2026")
banjir <- read_excel("data_banjir_pulausumatera.xlsx", sheet = "2010")
# Penamaan kolom agar konsisten dengan analisis selanjutnya
colnames(banjir) <- c("PROV", "TAHUN", "Y", "X1", "X2", "X3", "X4", "X5")
banjir$PROV <- as.factor(toupper(banjir$PROV))
banjir$TAHUN <- as.factor(banjir$TAHUN)
# Pre-Processing Data
# Cek Missing Value
print(colSums(is.na(banjir)))## PROV TAHUN Y X1 X2 X3 X4 X5
## 0 0 0 0 0 0 0 0# Cek Duplikasi Data
# Menghitung jumlah total baris yang duplikat
jumlah_duplikasi <- sum(duplicated(banjir))
print(paste("Jumlah data duplikat:", jumlah_duplikasi))## [1] "Jumlah data duplikat: 0"# Menampilkan baris data yang duplikat (jika ada)
if (jumlah_duplikasi > 0) {
print("Detail data yang terduplikasi:")
print(banjir[duplicated(banjir), ])
} else {
print("Tidak ditemukan data duplikat (bersih).")
}## [1] "Tidak ditemukan data duplikat (bersih)."3. Eksplorasi Data
Melakukan visualisasi untuk melihat pola keragaman kejadian banjir berdasarkan waktu dan antar provinsi, serta melihat korelasi antar variabel.
# 3. Eksplorasi Data ===========================================================
summary_stats <- banjir %>%
dplyr::select("Y", "X1", "X2", "X3", "X4","X5") %>%
summarise_all(list(
Min = ~min(.),
Q1 = ~quantile(., 0.25),
Median = ~median(.),
Mean = ~mean(.),
Q3 = ~quantile(., 0.75),
Max = ~max(.),
SD = ~sd(.)
)) %>%
pivot_longer(everything(),
names_to = c("Variable", "Statistic"),
names_sep = "_") %>%
pivot_wider(names_from = Statistic, values_from = value)
kable(summary_stats)| Variable | Min | Q1 | Median | Mean | Q3 | Max | SD |
|---|---|---|---|---|---|---|---|
| Y | 0.00 | 8.000 | 19.000 | 26.31333 | 36.0000 | 120.00 | 25.04094 |
| X1 | -9941.50 | 1826.150 | 5738.400 | 17635.24680 | 16856.1000 | 290777.00 | 38246.46864 |
| X2 | 116.33 | 5324.225 | 11004.000 | 24739.22720 | 22428.2825 | 225316.06 | 44814.03000 |
| X3 | 7.40 | 233.800 | 419.315 | 730.70507 | 1073.9500 | 3408.68 | 737.98993 |
| X4 | 91.50 | 191.160 | 231.090 | 233.33880 | 274.5475 | 435.67 | 64.77882 |
| X5 | 62.00 | 82.000 | 94.500 | 131.57333 | 194.5000 | 281.00 | 67.43273 |
# Keragaman Antar Waktu
ggplot(data = banjir, aes(x = TAHUN, y = Y)) +
geom_line() +
labs(x = "Tahun", y = "Jumlah Kejadian Banjir") +
theme(legend.position = "none") +
theme_bw()
# Keragaman Antar Individu
gplots::plotmeans(Y ~ `PROV`, main="Keragaman Antar Kejadian Banjir",
banjir, n.label = F)## Warning in arrows(x, li, x, pmax(y - gap, li), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped## Warning in arrows(x, ui, x, pmin(y + gap, ui), col = barcol, lwd = lwd, :
## zero-length arrow is of indeterminate angle and so skipped
# Di sini saya asumsikan continuous agar sumbu X rapi
banjir$TAHUN <- as.numeric(as.character(banjir$TAHUN))
# Plotting
p <- ggplot(data = banjir, aes(x = TAHUN, y = Y)) +
# 1. DATA MENTAH (BACKGROUND)
geom_point(position = position_jitter(width = 0.2),
alpha = 0.3,
color = "grey60",
size = 1.5) +
# 2. GARIS RATA-RATA (CENTRAL TENDENCY)
stat_summary(fun = mean,
geom = "line",
color = "black",
size = 0.8) +
# 3. TITIK RATA-RATA
stat_summary(fun = mean,
geom = "point",
shape = 18, # Shape diamond (formal)
size = 3,
color = "black") +
# 4. ERROR BARS (VARIABILITAS) - Sangat disukai reviewer
stat_summary(fun.data = mean_se,
geom = "errorbar",
width = 0.2,
size = 0.6) +
# 5. PEMFORMATAN LABEL & SKALA
scale_x_continuous(breaks = unique(banjir$TAHUN)) +
labs(x = "Tahun Pengamatan",
y = "Frekuensi Kejadian Banjir",
caption = "Sumber: [Sebutkan Data Anda Disini]\nError bars menunjukkan Standard Error of Mean (SEM)") +
# 6. TEMA ILMIAH (PUBLICATION READY THEME)
theme_bw() + # Dasar hitam putih
theme(
text = element_text(family = "serif"),
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10, color = "black"),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(), # Hilangkan minor grid
axis.line = element_line(colour = "black"), # Pertegas garis sumbu
plot.caption = element_text(hjust = 0, face = "italic", size = 9)
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(p)
# Correlation Plot
corrplot(cor(banjir[-c(1:2)]), method="color", type = "upper", tl.pos = 'tp')
corrplot(cor(banjir[-c(1:2)]), method="number", diag=F, add=T, type = "lower",
tl.pos = 'n', cl.pos = 'n')
4. Uji Multikolinearitas
Pengecekan masalah multikolinearitas antar variabel bebas menggunakan nilai VIF (Variance Inflation Factor).
# 4. Seleksi Peubah Multikolinearitas ==========================================
check_model <- lm(Y ~ X1 + X2 + X3 + X4 + X5, banjir)
data.frame(t(vif(check_model)))5. Pemodelan Data Panel (CEM & FEM)
Estimasi parameter menggunakan pendekatan Common Effect Model (CEM) dan Fixed Effect Model (FEM) dengan berbagai efek (Individu, Waktu, Two-ways).
# 5. Common-Effect Model =======================================================
# Model CEM
cem <- plm(Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, model = "pooling")
summary(cem)## Pooling Model
##
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, model = "pooling")
##
## Balanced Panel: n = 10, T = 15, N = 150
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -42.1958 -15.7700 -4.2596 7.6629 85.5516
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 7.8701e+00 9.3212e+00 0.8443 0.3998899
## X1 -5.3187e-05 5.8736e-05 -0.9055 0.3667030
## X2 -1.3782e-04 4.6543e-05 -2.9611 0.0035862 **
## X3 1.1399e-02 2.9160e-03 3.9089 0.0001422 ***
## X4 5.3063e-02 3.0217e-02 1.7561 0.0812027 .
## X5 1.5809e-02 3.1143e-02 0.5076 0.6124811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 93430
## Residual Sum of Squares: 78850
## R-Squared: 0.15605
## Adj. R-Squared: 0.12675
## F-statistic: 5.32541 on 5 and 144 DF, p-value: 0.00016048# 6. Fixed-Effect Model ========================================================
# Model FEM Individual
fem_ind <- plm(Y ~ X1 + X2 + X3 + X4 + X5, index = c("PROV", "TAHUN"),
model = "within", effect= "individual", data = banjir)
summary(fem_ind)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, effect = "individual",
## model = "within", index = c("PROV", "TAHUN"))
##
## Balanced Panel: n = 10, T = 15, N = 150
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -33.67095 -8.84941 0.11677 7.40798 57.77951
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## X1 -1.8483e-05 4.0638e-05 -0.4548 0.6499626
## X2 6.4756e-05 4.5498e-05 1.4233 0.1569622
## X3 3.0755e-02 8.0666e-03 3.8127 0.0002082 ***
## X4 1.1405e-01 2.7293e-02 4.1786 5.237e-05 ***
## X5 5.9241e-01 1.3487e-01 4.3925 2.246e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 45751
## Residual Sum of Squares: 30330
## R-Squared: 0.33706
## Adj. R-Squared: 0.26831
## F-statistic: 13.7275 on 5 and 135 DF, p-value: 7.7981e-11summary(fixef(fem_ind, effect="individual"))## Estimate Std. Error t-value Pr(>|t|)
## ACEH -25.533 13.960 -1.8290 0.0696 .
## BENGKULU -87.951 15.413 -5.7062 7.005e-08 ***
## JAMBI -71.540 13.824 -5.1752 8.060e-07 ***
## KEPULAUAN BANGKA BELITUNG -78.137 13.470 -5.8008 4.468e-08 ***
## KEPULAUAN RIAU -173.928 35.015 -4.9672 2.018e-06 ***
## LAMPUNG -152.341 33.351 -4.5678 1.099e-05 ***
## RIAU -126.550 26.082 -4.8520 3.321e-06 ***
## SUMATERA BARAT -88.869 19.713 -4.5081 1.405e-05 ***
## SUMATERA SELATAN -78.434 16.740 -4.6855 6.734e-06 ***
## SUMATERA UTARA -136.649 29.708 -4.5997 9.634e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model FEM Individual
fem_time <- plm(Y ~ X1 + X2 + X3 + X4 + X5, index = c("PROV", "TAHUN"),
model = "within", effect= "time", data = banjir)
summary(fem_time)## Oneway (time) effect Within Model
##
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, effect = "time",
## model = "within", index = c("PROV", "TAHUN"))
##
## Balanced Panel: n = 10, T = 15, N = 150
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -37.0485 -14.0108 -2.8886 8.0758 77.7908
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## X1 -1.4771e-05 6.1834e-05 -0.2389 0.811573
## X2 -1.0326e-04 5.0451e-05 -2.0467 0.042702 *
## X3 8.5406e-03 3.0264e-03 2.8221 0.005522 **
## X4 3.1019e-02 3.2625e-02 0.9508 0.343483
## X5 -2.7843e-03 3.1279e-02 -0.0890 0.929209
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 74769
## Residual Sum of Squares: 67998
## R-Squared: 0.090566
## Adj. R-Squared: -0.042352
## F-statistic: 2.5892 on 5 and 130 DF, p-value: 0.02875summary(fixef(fem_time, effect="time"))## Estimate Std. Error t-value Pr(>|t|)
## 2010 13.1029 12.8241 1.0217 0.308803
## 2011 4.9946 11.1656 0.4473 0.655389
## 2012 15.8313 11.4312 1.3849 0.168450
## 2013 9.8938 12.2266 0.8092 0.419877
## 2014 10.2313 11.0018 0.9300 0.354113
## 2015 9.8999 11.6416 0.8504 0.396675
## 2016 9.3849 12.0641 0.7779 0.438028
## 2017 7.8554 12.5719 0.6248 0.533174
## 2018 14.6589 12.2677 1.1949 0.234295
## 2019 9.8734 11.8697 0.8318 0.407039
## 2020 14.9420 12.3343 1.2114 0.227935
## 2021 32.8561 12.8340 2.5601 0.011608 *
## 2022 22.8874 13.1786 1.7367 0.084806 .
## 2023 33.8614 12.0868 2.8015 0.005864 **
## 2024 29.9688 13.1550 2.2781 0.024351 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model FEM Individual
fem_twoways <- plm(Y ~ X1 + X2 + X3 + X4 + X5, index = c("PROV", "TAHUN"),
model = "within", effect= "twoways", data = banjir)
summary(fem_twoways )## Twoways effects Within Model
##
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, effect = "twoways",
## model = "within", index = c("PROV", "TAHUN"))
##
## Balanced Panel: n = 10, T = 15, N = 150
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -30.80961 -8.04997 -0.14031 8.21278 56.02914
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## X1 -1.0900e-05 4.1367e-05 -0.2635 0.792621
## X2 8.3139e-06 4.8060e-05 0.1730 0.862948
## X3 1.1775e-02 1.0016e-02 1.1756 0.242050
## X4 1.0268e-01 3.2307e-02 3.1783 0.001881 **
## X5 4.4507e-02 2.4062e-01 0.1850 0.853562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 27090
## Residual Sum of Squares: 24399
## R-Squared: 0.099334
## Adj. R-Squared: -0.10908
## F-statistic: 2.669 on 5 and 121 DF, p-value: 0.025186summary(fixef(fem_twoways , effect="twoways"))## Estimate
## ACEH-2010 26.3305700
## ACEH-2011 23.1910797
## ACEH-2012 31.6055818
## ACEH-2013 26.9148115
## ACEH-2014 29.1415433
## ACEH-2015 32.0064124
## ACEH-2016 26.6575150
## ACEH-2017 25.4243314
## ACEH-2018 33.5629996
## ACEH-2019 30.8098126
## ACEH-2020 33.7741216
## ACEH-2021 50.3042083
## ACEH-2022 39.0737233
## ACEH-2023 53.4086043
## ACEH-2024 46.1518057
## BENGKULU-2010 -32.9855767
## BENGKULU-2011 -36.1250669
## BENGKULU-2012 -27.7105648
## BENGKULU-2013 -32.4013351
## BENGKULU-2014 -30.1746033
## BENGKULU-2015 -27.3097343
## BENGKULU-2016 -32.6586317
## BENGKULU-2017 -33.8918152
## BENGKULU-2018 -25.7531471
## BENGKULU-2019 -28.5063341
## BENGKULU-2020 -25.5420250
## BENGKULU-2021 -9.0119384
## BENGKULU-2022 -20.2424234
## BENGKULU-2023 -5.9075424
## BENGKULU-2024 -13.1643410
## JAMBI-2010 -21.7214970
## JAMBI-2011 -24.8609872
## JAMBI-2012 -16.4464851
## JAMBI-2013 -21.1372554
## JAMBI-2014 -18.9105236
## JAMBI-2015 -16.0456546
## JAMBI-2016 -21.3945520
## JAMBI-2017 -22.6277355
## JAMBI-2018 -14.4890673
## JAMBI-2019 -17.2422544
## JAMBI-2020 -14.2779453
## JAMBI-2021 2.2521414
## JAMBI-2022 -8.9783437
## JAMBI-2023 5.3565374
## JAMBI-2024 -1.9002613
## KEPULAUAN BANGKA BELITUNG-2010 -31.9175102
## KEPULAUAN BANGKA BELITUNG-2011 -35.0570004
## KEPULAUAN BANGKA BELITUNG-2012 -26.6424983
## KEPULAUAN BANGKA BELITUNG-2013 -31.3332686
## KEPULAUAN BANGKA BELITUNG-2014 -29.1065368
## KEPULAUAN BANGKA BELITUNG-2015 -26.2416678
## KEPULAUAN BANGKA BELITUNG-2016 -31.5905652
## KEPULAUAN BANGKA BELITUNG-2017 -32.8237487
## KEPULAUAN BANGKA BELITUNG-2018 -24.6850805
## KEPULAUAN BANGKA BELITUNG-2019 -27.4382676
## KEPULAUAN BANGKA BELITUNG-2020 -24.4739585
## KEPULAUAN BANGKA BELITUNG-2021 -7.9438718
## KEPULAUAN BANGKA BELITUNG-2022 -19.1743569
## KEPULAUAN BANGKA BELITUNG-2023 -4.8394758
## KEPULAUAN BANGKA BELITUNG-2024 -12.0962745
## KEPULAUAN RIAU-2010 -40.1239724
## KEPULAUAN RIAU-2011 -43.2634627
## KEPULAUAN RIAU-2012 -34.8489606
## KEPULAUAN RIAU-2013 -39.5397309
## KEPULAUAN RIAU-2014 -37.3129991
## KEPULAUAN RIAU-2015 -34.4481300
## KEPULAUAN RIAU-2016 -39.7970274
## KEPULAUAN RIAU-2017 -41.0302110
## KEPULAUAN RIAU-2018 -32.8915428
## KEPULAUAN RIAU-2019 -35.6447298
## KEPULAUAN RIAU-2020 -32.6804208
## KEPULAUAN RIAU-2021 -16.1503341
## KEPULAUAN RIAU-2022 -27.3808191
## KEPULAUAN RIAU-2023 -13.0459381
## KEPULAUAN RIAU-2024 -20.3027367
## LAMPUNG-2010 -20.0714366
## LAMPUNG-2011 -23.2109269
## LAMPUNG-2012 -14.7964248
## LAMPUNG-2013 -19.4871951
## LAMPUNG-2014 -17.2604633
## LAMPUNG-2015 -14.3955942
## LAMPUNG-2016 -19.7444916
## LAMPUNG-2017 -20.9776751
## LAMPUNG-2018 -12.8390070
## LAMPUNG-2019 -15.5921940
## LAMPUNG-2020 -12.6278849
## LAMPUNG-2021 3.9022017
## LAMPUNG-2022 -7.3282833
## LAMPUNG-2023 7.0065977
## LAMPUNG-2024 -0.2502009
## RIAU-2010 -40.0599813
## RIAU-2011 -43.1994716
## RIAU-2012 -34.7849695
## RIAU-2013 -39.4757398
## RIAU-2014 -37.2490080
## RIAU-2015 -34.3841389
## RIAU-2016 -39.7330363
## RIAU-2017 -40.9662199
## RIAU-2018 -32.8275517
## RIAU-2019 -35.5807387
## RIAU-2020 -32.6164297
## RIAU-2021 -16.0863430
## RIAU-2022 -27.3168280
## RIAU-2023 -12.9819470
## RIAU-2024 -20.2387456
## SUMATERA BARAT-2010 -14.9192267
## SUMATERA BARAT-2011 -18.0587169
## SUMATERA BARAT-2012 -9.6442148
## SUMATERA BARAT-2013 -14.3349851
## SUMATERA BARAT-2014 -12.1082533
## SUMATERA BARAT-2015 -9.2433843
## SUMATERA BARAT-2016 -14.5922817
## SUMATERA BARAT-2017 -15.8254652
## SUMATERA BARAT-2018 -7.6867970
## SUMATERA BARAT-2019 -10.4399841
## SUMATERA BARAT-2020 -7.4756750
## SUMATERA BARAT-2021 9.0544117
## SUMATERA BARAT-2022 -2.1760734
## SUMATERA BARAT-2023 12.1588077
## SUMATERA BARAT-2024 4.9020090
## SUMATERA SELATAN-2010 -13.5026882
## SUMATERA SELATAN-2011 -16.6421784
## SUMATERA SELATAN-2012 -8.2276763
## SUMATERA SELATAN-2013 -12.9184466
## SUMATERA SELATAN-2014 -10.6917148
## SUMATERA SELATAN-2015 -7.8268458
## SUMATERA SELATAN-2016 -13.1757432
## SUMATERA SELATAN-2017 -14.4089267
## SUMATERA SELATAN-2018 -6.2702585
## SUMATERA SELATAN-2019 -9.0234456
## SUMATERA SELATAN-2020 -6.0591365
## SUMATERA SELATAN-2021 10.4709502
## SUMATERA SELATAN-2022 -0.7595349
## SUMATERA SELATAN-2023 13.5753462
## SUMATERA SELATAN-2024 6.3185475
## SUMATERA UTARA-2010 -7.8306517
## SUMATERA UTARA-2011 -10.9701420
## SUMATERA UTARA-2012 -2.5556399
## SUMATERA UTARA-2013 -7.2464102
## SUMATERA UTARA-2014 -5.0196784
## SUMATERA UTARA-2015 -2.1548093
## SUMATERA UTARA-2016 -7.5037067
## SUMATERA UTARA-2017 -8.7368903
## SUMATERA UTARA-2018 -0.5982221
## SUMATERA UTARA-2019 -3.3514091
## SUMATERA UTARA-2020 -0.3871001
## SUMATERA UTARA-2021 16.1429866
## SUMATERA UTARA-2022 4.9125016
## SUMATERA UTARA-2023 19.2473826
## SUMATERA UTARA-2024 11.99058396. Evaluasi Model FEM
Melakukan uji spesifikasi untuk menentukan efek yang signifikan pada FEM serta membandingkan performa model berdasarkan kriteria kebaikan model (SSE, AIC, MAPE, R-Squared).
# 7. Uji Signifikasi Pengaruh FEM ==============================================
# Uji pengaruh individu
plmtest(fem_twoways, type = "bp", effect = "individual")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 238.01, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsplmtest(fem_twoways, type = "bp", effect = "time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 0.49453, df = 1, p-value = 0.4819
## alternative hypothesis: significant effectsplmtest(fem_twoways, type = "bp", effect = "twoways")##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 238.5, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects# Nilai Kebaikan Model
# Sum Squared Error
dsse <- data.frame(Individu=sum(fem_ind$residuals^2),
Time=sum(fem_time$residuals^2),
Twoways=sum(fem_twoways$residuals^2))
# AIC
lsdv_ind <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + PROV, banjir)
lsdv_time <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + TAHUN, banjir)
lsdv_twoways <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + PROV + TAHUN, banjir)
daic <- data.frame(Individu=AIC(lsdv_ind),Time=AIC(lsdv_time),
Twoways=AIC(lsdv_twoways))
# MAPE
mape <- function(actual, forecast) {
mean(abs((actual - forecast) / actual)) * 100
}
dmape <- data.frame(Individu=mape(banjir$Y,predict(fem_ind)),
Time=mape(banjir$Y,predict(fem_time)),
Twoways=mape(banjir$Y,predict(fem_twoways)))
# BIC
dbic <- data.frame(Individu=BIC(lsdv_ind),Time=BIC(lsdv_time),
Twoways=BIC(lsdv_twoways))
# R-squared
ind <- summary(fem_ind)
time <- summary(fem_time)
tways <- summary(fem_twoways)
drsq <- data.frame(Individu=ind$r.squared,Time=time$r.squared,
Twoways=tways$r.squared)
# Perbandingan
compare <- t(rbind(dsse,daic,dbic,dmape,drsq))
colnames(compare) <- c("SSE","AIC","BIC", "MAPE","R-Squared","Adj R-Squared")
compare## SSE AIC BIC MAPE R-Squared Adj R-Squared
## Individu 30329.96 1254.070 1302.240 Inf 0.33705789 0.26830834
## Time 67997.78 1370.530 1394.615 Inf 0.09056555 -0.04235179
## Twoways 24398.65 1247.023 1298.204 Inf 0.09933384 -0.10908478# FEM vs CEM
pooltest(cem, fem_ind)##
## F statistic
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## F = 23.996, df1 = 9, df2 = 135, p-value < 2.2e-16
## alternative hypothesis: unstability7. Random Effect Model (REM)
Estimasi model menggunakan pendekatan Random Effect dan melakukan uji Lagrange Multiplier serta uji Hausman untuk membandingkan FEM dan REM.
# 8. Random-Effect Model =======================================================
# Pendugaan Parameter
rem_gls <- plm(Y ~ X1 + X2 + X3 + X4 + X5, data = banjir,
index = c("PROV", "TAHUN"),
effect = "individual", model = "random",
random.method = "nerlove")
summary(rem_gls)## Oneway (individual) effect Random Effect Model
## (Nerlove's transformation)
##
## Call:
## plm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = banjir, effect = "individual",
## model = "random", random.method = "nerlove", index = c("PROV",
## "TAHUN"))
##
## Balanced Panel: n = 10, T = 15, N = 150
##
## Effects:
## var std.dev share
## idiosyncratic 202.20 14.22 0.093
## individual 1975.64 44.45 0.907
## theta: 0.9177
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -35.3964 -8.5495 -1.8589 7.3066 58.2127
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -8.2025e+01 2.3087e+01 -3.5529 0.000381 ***
## X1 -1.9621e-05 3.9927e-05 -0.4914 0.623126
## X2 3.7422e-05 4.3442e-05 0.8614 0.389008
## X3 2.9209e-02 7.4117e-03 3.9409 8.117e-05 ***
## X4 1.1617e-01 2.6935e-02 4.3129 1.611e-05 ***
## X5 4.5077e-01 1.1459e-01 3.9337 8.363e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 46074
## Residual Sum of Squares: 31796
## R-Squared: 0.3099
## Adj. R-Squared: 0.28594
## Chisq: 64.6647 on 5 DF, p-value: 1.3153e-12# Perbandingan Kebaikan REM
# Efek individu
plmtest(rem_gls,type = "bp", effect="individu")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 238.01, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects# Efek time
plmtest(rem_gls,type = "bp", effect="time")##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 0.49453, df = 1, p-value = 0.4819
## alternative hypothesis: significant effects# Efek twoways
plmtest(rem_gls,type = "bp", effect="twoways")##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 238.5, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effects# FEM vs REM
#fem dengan rem gls
phtest(fem_ind, rem_gls)##
## Hausman Test
##
## data: Y ~ X1 + X2 + X3 + X4 + X5
## chisq = 4.8307, df = 5, p-value = 0.4369
## alternative hypothesis: one model is inconsistent8. Pemilihan Model Terbaik
Estimasi model menggunakan pendekatan Random Effect dan melakukan uji Lagrange Multiplier serta uji Hausman untuk membandingkan FEM dan REM.
# =========================================================================
# UJI PEMILIHAN MODEL (CHOW, HAUSMAN, LM) & TABEL REKAPITULASI
# =========================================================================
# A. Uji Chow (Likelihood Ratio Test): Membandingkan CEM vs FEM
# H0: Model CEM lebih baik
# H1: Model FEM lebih baik
chow_test <- pooltest(cem, fem_ind)
# B. Uji Hausman: Membandingkan FEM vs REM
# H0: Model REM lebih baik (Konsisten & Efisien)
# H1: Model FEM lebih baik (Konsisten)
hausman_test <- phtest(fem_ind, rem_gls)
# C. Uji LM (Lagrange Multiplier - Breusch-Pagan): Membandingkan CEM vs REM
# H0: Variansi error individual = 0 (Model CEM lebih baik)
# H1: Variansi error individual != 0 (Model REM lebih baik)
lm_test <- plmtest(cem, type = "bp", effect = "individual")
# --- FUNGSI UNTUK MENENTUKAN KESIMPULAN OTOMATIS ---
get_conclusion <- function(p_val, test_type) {
alpha <- 0.05
if (test_type == "Chow") {
return(ifelse(p_val < alpha, "Model FEM lebih baik", "Model CEM lebih baik"))
} else if (test_type == "Hausman") {
return(ifelse(p_val < alpha, "Model FEM lebih baik", "Model REM lebih baik"))
} else if (test_type == "LM") {
return(ifelse(p_val < alpha, "Model REM lebih baik", "Model CEM lebih baik"))
}
}
# --- MEMBUAT DATA FRAME HASIL ---
hasil_uji <- data.frame(
Uji = c("Chow (F-Test)", "Hausman", "Lagrange Multiplier (LM)"),
Statistik = c(
as.numeric(chow_test$statistic),
as.numeric(hausman_test$statistic),
as.numeric(lm_test$statistic)
),
P_Value = c(
as.numeric(chow_test$p.value),
as.numeric(hausman_test$p.value),
as.numeric(lm_test$p.value)
),
Kesimpulan = c(
get_conclusion(as.numeric(chow_test$p.value), "Chow"),
get_conclusion(as.numeric(hausman_test$p.value), "Hausman"),
get_conclusion(as.numeric(lm_test$p.value), "LM")
)
)
# --- FORMATTING TAMPILAN AGAR RAPI (OPSIONAL) ---
hasil_uji$Statistik <- round(hasil_uji$Statistik, 4)
# Mengubah P-value yang sangat kecil menjadi notasi ilmiah atau < 0.05
hasil_uji$P_Value_Str <- ifelse(hasil_uji$P_Value < 0.0001,
format(hasil_uji$P_Value, scientific = TRUE, digits = 3),
round(hasil_uji$P_Value, 4))
# Menampilkan Tabel Akhir (Bersih)
library(knitr) # Gunakan library ini jika ingin tampilan tabel cantik di console
tabel_final <- hasil_uji[, c("Uji", "Statistik", "P_Value_Str", "Kesimpulan")]
colnames(tabel_final)[3] <- "P-value" # Ubah nama kolom biar sesuai contoh
print("=== TABEL RINGKASAN PEMILIHAN MODEL PANEL ===")## [1] "=== TABEL RINGKASAN PEMILIHAN MODEL PANEL ==="kable(tabel_final, align = "c")| Uji | Statistik | P-value | Kesimpulan |
|---|---|---|---|
| Chow (F-Test) | 23.9962 | 4.52e-24 | Model FEM lebih baik |
| Hausman | 4.8307 | 0.4369 | Model REM lebih baik |
| Lagrange Multiplier (LM) | 238.0099 | 1.07e-53 | Model REM lebih baik |
9. Uji Diagnostik Asumsi & Spasial
Pemeriksaan asumsi klasik (Normalitas, Heteroskedastisitas, Autokorelasi) serta pengujian dependensi spasial (Cross-sectional Dependence) untuk melegitimasi penggunaan model spasial.
# 9. UJi Diagnostik Sisaan =====================================================
sum_rem <- summary(rem_gls)
# Menghitung VIF
vif_values <- tryCatch(vif(rem_gls), error = function(e) NULL)
if (!is.null(vif_values)) {
if (length(coef(rem_gls)) > length(vif_values)) {
vif_vec <- c(NA, vif_values)
} else {
vif_vec <- vif_values
}
} else {
vif_vec <- rep(NA, length(coef(rem_gls)))
}
# --- 2. Membuat Tabel Estimasi (Koefisien & VIF) ---
tabel_estimasi <- data.frame(
Variabel = rownames(sum_rem$coefficients),
Koefisien = round(sum_rem$coefficients[,1], 4),
Std_Error = round(sum_rem$coefficients[,2], 4),
P_Value = round(sum_rem$coefficients[,4], 4),
VIF = round(vif_vec, 2)
)
tabel_estimasi$Signifikansi <- ifelse(tabel_estimasi$P_Value < 0.05,
"Signifikan (5%)",
ifelse(tabel_estimasi$P_Value < 0.1,
"Signifikan (10%)", "Tidak Signifikan"))
tabel_estimasi$Status_VIF <- ifelse(is.na(tabel_estimasi$VIF), "-",
ifelse(tabel_estimasi$VIF < 10,
"Aman", "Multikolinearitas"))
# --- 3. Melakukan Uji Diagnostik & Membuat Tabel ---
# A. Uji Normalitas (Kolmogorov-Smirnov)
uji_norm <- ks.test(residuals(rem_gls), "pnorm",
mean=mean(residuals(rem_gls)),
sd=sd(residuals(rem_gls)))
# B. Uji Heteroskedastisitas (Breusch-Pagan)
uji_hetero <- bptest(rem_gls, studentize = TRUE)
# C. Uji Autokorelasi Serial (Waktu - Breusch-Godfrey)
uji_auto_serial <- pbgtest(rem_gls)
# D. UJI BARU: Uji Dependensi Spasial / Cross-sectional Dependence (Pesaran CD)
# H0: Tidak ada dependensi antar wilayah (Cross-sectional independence)
# H1: Ada dependensi antar wilayah (Indikasi Spasial)->P < 0.05 BAGUS untuk GWPR
uji_spasial <- pcdtest(rem_gls, test = "cd")
# Mengambil R-Squared
r_sq <- sum_rem$r.squared[1]
# Menyusun Tabel Diagnostik
tabel_diagnostik <- data.frame(
Uji_Asumsi = c("Normalitas (KS)",
"Heteroskedastisitas (BP)",
"Autokorelasi Serial (Waktu)",
"Dependensi Spasial (Pesaran CD)", # Baris Baru
"R-Squared Overall"),
Nilai_P = c(round(uji_norm$p.value, 4),
round(uji_hetero$p.value, 4),
round(uji_auto_serial$p.value, 4),
round(uji_spasial$p.value, 4), # Nilai P Spasial
NA),
Statistik = c(round(uji_norm$statistic, 4),
round(uji_hetero$statistic, 4),
round(uji_auto_serial$statistic, 4),
round(uji_spasial$statistic, 4), # Statistik Spasial
round(r_sq, 4)),
Keputusan = c(ifelse(uji_norm$p.value > 0.05,
"Normal (Terpenuhi)", "Tidak Normal"),
ifelse(uji_hetero$p.value > 0.05,
"Homoskedastis (Terpenuhi)",
"Heteroskedastis (Cek Spasial)"),
ifelse(uji_auto_serial$p.value > 0.05,
"Bebas Autokorelasi Waktu", "Ada Autokorelasi Waktu"),
ifelse(uji_spasial$p.value < 0.05,
"Ada Efek Spasial (Lanjut GWPR)",
"Tidak Ada Efek Spasial"),
paste0(round(r_sq * 100, 2), "%"))
)
# --- 4. Cetak Output dengan Kable ---
cat("\n=== RINGKASAN ESTIMASI RANDOM EFFECT MODEL (REM) ===\n")##
## === RINGKASAN ESTIMASI RANDOM EFFECT MODEL (REM) ===print(kable(tabel_estimasi, row.names=FALSE, align="c",
caption = "Koefisien dan Uji Multikolinearitas"))##
##
## Table: Koefisien dan Uji Multikolinearitas
##
## | Variabel | Koefisien | Std_Error | P_Value | VIF | Signifikansi | Status_VIF |
## |:-----------:|:---------:|:---------:|:-------:|:----:|:----------------:|:----------:|
## | (Intercept) | -82.0254 | 23.0869 | 0.0004 | NA | Signifikan (5%) | - |
## | X1 | 0.0000 | 0.0000 | 0.6231 | 1.20 | Tidak Signifikan | Aman |
## | X2 | 0.0000 | 0.0000 | 0.3890 | 1.57 | Tidak Signifikan | Aman |
## | X3 | 0.0292 | 0.0074 | 0.0001 | 1.35 | Signifikan (5%) | Aman |
## | X4 | 0.1162 | 0.0269 | 0.0000 | 1.05 | Signifikan (5%) | Aman |
## | X5 | 0.4508 | 0.1146 | 0.0001 | 1.24 | Signifikan (5%) | Aman |cat("\n=== DIAGNOSTIK ASUMSI KLASIK & SPASIAL ===\n")##
## === DIAGNOSTIK ASUMSI KLASIK & SPASIAL ===print(kable(tabel_diagnostik, row.names=FALSE, align="l",
caption = "Hasil Uji Diagnostik"))##
##
## Table: Hasil Uji Diagnostik
##
## |Uji_Asumsi |Nilai_P |Statistik |Keputusan |
## |:-------------------------------|:-------|:---------|:------------------------------|
## |Normalitas (KS) |0.3138 |0.0785 |Normal (Terpenuhi) |
## |Heteroskedastisitas (BP) |0.5460 |4.0239 |Homoskedastis (Terpenuhi) |
## |Autokorelasi Serial (Waktu) |0.0014 |36.6224 |Ada Autokorelasi Waktu |
## |Dependensi Spasial (Pesaran CD) |0.0204 |2.3181 |Ada Efek Spasial (Lanjut GWPR) |
## |R-Squared Overall |NA |0.3099 |30.99% |# Histogram
ggplot(as.data.frame(rem_gls$residuals), aes(x = rem_gls$residuals)) +
geom_histogram(aes(y = after_stat(density)), color = "white",
fill = "steelblue") +
geom_density(color = "red", linewidth = 1) +
theme_minimal()## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
10. Pemodelan GWPR
Memuat data koordinat wilayah, menentukan bandwidth optimal (Adaptive/Fixed) dengan berbagai fungsi kernel (Bisquare, Gaussian, Exponential), dan membangun model GWPR.
# 10 Geographically Weighted Panel Regression ==================================
# Input Latitude Longitude
wilayah <- fromJSON("https://raw.githubusercontent.com/yusufsyaifudin/wilayah-indonesia/master/data/list_of_area/provinces.json")
wilayah$name <- toupper(wilayah$name)
data.sp.GWPR <- banjir %>%
left_join(wilayah %>% select(name, latitude, longitude),
by = c("PROV" = "name")) %>%
rename(LAT = latitude, LONG = longitude) %>%
select(PROV, TAHUN, LAT, LONG, everything())
# Modeling GWPR
coordinates(data.sp.GWPR)=3:4
class(data.sp.GWPR)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"# Menentukan Fungsi Pembobot Spasial Terbaik
# Fungsi kernel: ((gaussian, eksponensial, bisquare, Tricube) (adaptive dan fixed))
# --- KELOMPOK ADAPTIVE ---
# 1. Adaptive Bisquare
bwd.GWPR.bisquare.ad <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "bisquare", adaptive=T)## Adaptive bandwidth: 100 CV score: 43636.5
## Adaptive bandwidth: 70 CV score: 41819.79
## Adaptive bandwidth: 50 CV score: 38156.24
## Adaptive bandwidth: 39 CV score: 36535.15
## Adaptive bandwidth: 31 CV score: 36535.15hasil.GWPR.bisquare.ad <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.bisquare.ad,
kernel = "bisquare", adaptive=T)
# 2. Adaptive Gaussian
bwd.GWPR.gaussian.ad <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "gaussian", adaptive=T)## Adaptive bandwidth: 100 CV score: 70008.75
## Adaptive bandwidth: 70 CV score: 65964.49
## Adaptive bandwidth: 50 CV score: 58937.44
## Adaptive bandwidth: 39 CV score: 55737.64
## Adaptive bandwidth: 31 CV score: 55737.64hasil.GWPR.gaussian.ad <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.gaussian.ad,
kernel = "gaussian", adaptive=T)
# 3. Adaptive Exponential
bwd.GWPR.exponential.ad <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "exponential", adaptive=T)## Adaptive bandwidth: 100 CV score: 62286.85
## Adaptive bandwidth: 70 CV score: 59273.64
## Adaptive bandwidth: 50 CV score: 55106.61
## Adaptive bandwidth: 39 CV score: 53314.57
## Adaptive bandwidth: 31 CV score: 53314.57hasil.GWPR.exponential.ad <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.exponential.ad,
kernel = "exponential", adaptive=T)
# 4. Adaptive Tricube (BARU)
bwd.GWPR.tricube.ad <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "tricube", adaptive=T)## Adaptive bandwidth: 100 CV score: 44122.37
## Adaptive bandwidth: 70 CV score: 41745.24
## Adaptive bandwidth: 50 CV score: 37500.37
## Adaptive bandwidth: 39 CV score: 36631.9
## Adaptive bandwidth: 31 CV score: 36631.9hasil.GWPR.tricube.ad <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.tricube.ad,
kernel = "tricube", adaptive=T)
# --- KELOMPOK FIXED ---
# 5. Bisquare
bwd.GWPR.bisquare.fix <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "bisquare", adaptive=F)## Fixed bandwidth: 7.873791 CV score: 48073
## Fixed bandwidth: 4.867243 CV score: 39069.78
## Fixed bandwidth: 3.009095 CV score: 34268.34
## Fixed bandwidth: 1.860696 CV score: 26792.34
## Fixed bandwidth: 1.150947 CV score: 28458.12
## Fixed bandwidth: 2.299345 CV score: 24711.74
## Fixed bandwidth: 2.570446 CV score: 22971.01
## Fixed bandwidth: 2.737995 CV score: 23261.04
## Fixed bandwidth: 2.466895 CV score: 23288.38
## Fixed bandwidth: 2.634444 CV score: 23020.68
## Fixed bandwidth: 2.530893 CV score: 23055.61
## Fixed bandwidth: 2.594891 CV score: 22968.27
## Fixed bandwidth: 2.609999 CV score: 22981.67
## Fixed bandwidth: 2.585554 CV score: 22965.36
## Fixed bandwidth: 2.579783 CV score: 22965.93
## Fixed bandwidth: 2.58912 CV score: 22965.94
## Fixed bandwidth: 2.583349 CV score: 22965.36
## Fixed bandwidth: 2.581987 CV score: 22965.49
## Fixed bandwidth: 2.584191 CV score: 22965.33
## Fixed bandwidth: 2.584712 CV score: 22965.33
## Fixed bandwidth: 2.58387 CV score: 22965.33
## Fixed bandwidth: 2.58439 CV score: 22965.33
## Fixed bandwidth: 2.584513 CV score: 22965.33
## Fixed bandwidth: 2.584314 CV score: 22965.33
## Fixed bandwidth: 2.584437 CV score: 22965.33hasil.GWPR.bisquare.fix <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.bisquare.fix,
kernel = "bisquare", adaptive=F)
# 6. Fixed Gaussian
bwd.GWPR.gaussian.fix <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "gaussian", adaptive=F)## Fixed bandwidth: 7.873791 CV score: 69333.52
## Fixed bandwidth: 4.867243 CV score: 57038.7
## Fixed bandwidth: 3.009095 CV score: 46845.79
## Fixed bandwidth: 1.860696 CV score: 39337.6
## Fixed bandwidth: 1.150947 CV score: 34167.23
## Fixed bandwidth: 0.7122972 CV score: 24725.82
## Fixed bandwidth: 0.441197 CV score: 27344.92
## Fixed bandwidth: 0.8798464 CV score: 31464.35
## Fixed bandwidth: 0.6087462 CV score: 24466.57
## Fixed bandwidth: 0.5447481 CV score: 25102.2
## Fixed bandwidth: 0.6482991 CV score: 24144.46
## Fixed bandwidth: 0.6727442 CV score: 24113.23
## Fixed bandwidth: 0.6878521 CV score: 24229.29
## Fixed bandwidth: 0.663407 CV score: 24099.26
## Fixed bandwidth: 0.6576363 CV score: 24108.04
## Fixed bandwidth: 0.6669735 CV score: 24100.14
## Fixed bandwidth: 0.6612028 CV score: 24101.21
## Fixed bandwidth: 0.6647693 CV score: 24098.99
## Fixed bandwidth: 0.6656113 CV score: 24099.19
## Fixed bandwidth: 0.664249 CV score: 24099.01
## Fixed bandwidth: 0.6650909 CV score: 24099.03
## Fixed bandwidth: 0.6645706 CV score: 24098.98
## Fixed bandwidth: 0.6644477 CV score: 24098.99
## Fixed bandwidth: 0.6646465 CV score: 24098.98
## Fixed bandwidth: 0.6645236 CV score: 24098.98
## Fixed bandwidth: 0.6645996 CV score: 24098.98hasil.GWPR.gaussian.fix <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.gaussian.fix,
kernel = "gaussian", adaptive=F)
# 7. Fixed Exponential
bwd.GWPR.exponential.fix <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "exponential", adaptive=F)## Fixed bandwidth: 7.873791 CV score: 62204.64
## Fixed bandwidth: 4.867243 CV score: 54130.08
## Fixed bandwidth: 3.009095 CV score: 47025.43
## Fixed bandwidth: 1.860696 CV score: 41651.48
## Fixed bandwidth: 1.150947 CV score: 37707.6
## Fixed bandwidth: 0.7122972 CV score: 33565.85
## Fixed bandwidth: 0.441197 CV score: 28431.57
## Fixed bandwidth: 0.2736479 CV score: 25288.98
## Fixed bandwidth: 0.1700968 CV score: 28185.32
## Fixed bandwidth: 0.337646 CV score: 25118.05
## Fixed bandwidth: 0.377199 CV score: 25936.37
## Fixed bandwidth: 0.3132009 CV score: 25010.59
## Fixed bandwidth: 0.298093 CV score: 25063.39
## Fixed bandwidth: 0.3225381 CV score: 25020.24
## Fixed bandwidth: 0.3074302 CV score: 25021.55
## Fixed bandwidth: 0.3167674 CV score: 25010.09
## Fixed bandwidth: 0.3189716 CV score: 25012.33
## Fixed bandwidth: 0.3154051 CV score: 25009.69
## Fixed bandwidth: 0.3145632 CV score: 25009.81
## Fixed bandwidth: 0.3159254 CV score: 25009.76
## Fixed bandwidth: 0.3150835 CV score: 25009.71
## Fixed bandwidth: 0.3156038 CV score: 25009.71
## Fixed bandwidth: 0.3152823 CV score: 25009.69hasil.GWPR.exponential.fix <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.exponential.fix,
kernel = "exponential", adaptive=F)
# 8. Fixed Tricube (BARU)
bwd.GWPR.tricube.fix <- bw.gwr(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR, approach = "CV",
kernel = "tricube", adaptive=F)## Fixed bandwidth: 7.873791 CV score: 48442.01
## Fixed bandwidth: 4.867243 CV score: 39172.72
## Fixed bandwidth: 3.009095 CV score: 33207.67
## Fixed bandwidth: 1.860696 CV score: 26894.05
## Fixed bandwidth: 1.150947 CV score: 28458.12
## Fixed bandwidth: 2.299345 CV score: 24840.73
## Fixed bandwidth: 2.570446 CV score: 23057.9
## Fixed bandwidth: 2.737995 CV score: 23004.35
## Fixed bandwidth: 2.841546 CV score: 24976.5
## Fixed bandwidth: 2.673997 CV score: 22903.65
## Fixed bandwidth: 2.634444 CV score: 22912.8
## Fixed bandwidth: 2.698442 CV score: 22929.27
## Fixed bandwidth: 2.658889 CV score: 22899.05
## Fixed bandwidth: 2.649552 CV score: 22901.11
## Fixed bandwidth: 2.664659 CV score: 22899.69
## Fixed bandwidth: 2.655322 CV score: 22899.38
## Fixed bandwidth: 2.661093 CV score: 22899.13
## Fixed bandwidth: 2.657526 CV score: 22899.11
## Fixed bandwidth: 2.659731 CV score: 22899.06
## Fixed bandwidth: 2.658368 CV score: 22899.07
## Fixed bandwidth: 2.65921 CV score: 22899.05
## Fixed bandwidth: 2.659409 CV score: 22899.05
## Fixed bandwidth: 2.659088 CV score: 22899.05
## Fixed bandwidth: 2.659286 CV score: 22899.05
## Fixed bandwidth: 2.659163 CV score: 22899.05hasil.GWPR.tricube.fix <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
data = data.sp.GWPR,
bw = bwd.GWPR.tricube.fix,
kernel = "tricube", adaptive=F)11. Pemilihan Model Terbaik & Parameter Lokal
Membandingkan kebaikan model berdasarkan AIC/AICc, R-Squared, dan RSS, serta menampilkan parameter lokal dan signifikansinya per provinsi.
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.tricube.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,
hasil.GWPR.tricube.fix$GW.diagnostic$AIC)
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.tricube.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,
hasil.GWPR.tricube.fix$GW.diagnostic$gw.R2)
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.tricube.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,
hasil.GWPR.tricube.fix$GW.diagnostic$AICc)
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.tricube.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,
hasil.GWPR.tricube.fix$GW.diagnostic$gwR2.adj)
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.tricube.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,
hasil.GWPR.tricube.fix$GW.diagnostic$RSS.gw)
tabel <- cbind(RSS,R2,R2adj,AIC,AICc)
rownames(tabel) <- c("Adaptive Bisquare", "Adaptive Gaussian",
"Adaptive Exponential", "Adaptive Tricube", "Fixed Bisquare",
"Fixed Gaussian", "Fixed Exponential", "Fixed Tricube")
data.frame(tabel)# Pendugaan Parameter
parameter.GWPR <- as.data.frame(
hasil.GWPR.tricube.fix$SDF[1:10,1:6])[-c(7,8)]
parameter.GWPR %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, everything())parameter_R2 <- as.data.frame(
hasil.GWPR.tricube.fix$SDF[1:10,c(1:6,24)])[-8:-19]
# R2 setiap provinsi
parameter_R2 %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, Local_R2) %>%
arrange(-Local_R2)# 11 Signifikasi parameter (p-value) ===========================================
p_value <- gwr.t.adjust(hasil.GWPR.tricube.fix)$results$p
p_value <- as.data.frame(ifelse(p_value<=0.05, "Signifikan", "Tidak"))[1:10,]
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, everything())#Filter Provinsi yang mempunyai signifikansi pada peubah X1
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, X1_p) %>%
filter("Signifikan" == p_value$X1)#Filter Provinsi yang mempunyai signifikansi pada peubah X2
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, X2_p) %>%
filter("Signifikan" == p_value$X2)#Filter Provinsi yang mempunyai signifikansi pada peubah X3
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, X3_p) %>%
filter("Signifikan" == p_value$X3)#Filter Provinsi yang mempunyai signifikansi pada peubah X4
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, X4_p) %>%
filter("Signifikan" == p_value$X4)# Filter Provinsi yang mempunyai signifikansi pada peubah X5
p_value %>%
mutate(PROV=banjir$PROV[1:10]) %>%
select(PROV, X5_p) %>%
filter("Signifikan" == p_value$X5)# Lihat koefisien untuk provinsi tertentu (misal baris ke-1 adalah Sumatera Utara)
hasil.GWPR.tricube.fix$SDF[1,] ## coordinates Intercept X1 X2 X3 X4
## 1 (4.36855, 97.0253) 132.0107 -0.0005626677 -0.003256374 0.2081083 0.1654052
## X5 y yhat residual CV_Score Stud_residual Intercept_SE
## 1 -1.533336 36 35.21741 0.7825929 0 0.1102028 84.45181
## X1_SE X2_SE X3_SE X4_SE X5_SE Intercept_TV
## 1 0.0004214656 0.0009942214 0.1355392 0.07654525 1.120265 1.563148
## X1_TV X2_TV X3_TV X4_TV X5_TV Local_R2
## 1 -1.335027 -3.2753 1.53541 2.160882 -1.368726 0.9256041# Atau jika ingin melihat signifikansi (t-hitung) untuk memberi tanda bintang (*)
gwr.t.adjust(hasil.GWPR.tricube.fix)$results$t[1,]## Intercept_t X1_t X2_t X3_t X4_t X5_t
## 1.563148 -1.335027 -3.275300 1.535410 2.160882 -1.36872612. Uji Kesesuaian Model (Global vs Lokal)
Melakukan uji F untuk membandingkan apakah model GWPR (lokal) secara signifikan lebih baik daripada model global (REM).
# 12. Uji Kesesuaian Model ====================================================
# Uji F
rss_rem <- sum(rem_gls$residuals^2)
rss_gwpr <- hasil.GWPR.tricube.fix$GW.diagnostic$RSS.gw
df_rem <- rem_gls$df.residual
df_gwpr <- hasil.GWPR.tricube.fix$GW.diagnostic$edf
Fhit=(rss_gwpr/df_gwpr)/(rss_rem/df_rem)
Fhit## [1] 0.5599753pf(Fhit,df_gwpr,df_rem)## [1] 0.0008418435comparison <- data.frame(R2=c(hasil.GWPR.tricube.fix$GW.diagnostic$gw.R2,
summary(rem_gls)$r.squared[1]),
R2adj=c(hasil.GWPR.tricube.fix$GW.diagnostic$gwR2.adj,
summary(rem_gls)$r.squared[2]),
RSS=c(hasil.GWPR.tricube.fix$GW.diagnostic$RSS.gw,
sum(summary(rem_gls)$residuals^2)))
rownames(comparison) <- c("GWPR", "REM Two Ways")
comparison# ==============================================================================
# UJI ANOVA (UJI F) PERBANDINGAN MODEL GLOBAL VS GWPR
# ==============================================================================
# 1. Masukkan Nilai RSS dan DF (Gunakan Data dari Tabel 6 Anda agar Konsisten)
# ------------------------------------------------------------------------------
# Nilai ini diambil dari Tabel 6 yang Anda tunjukkan sebelumnya (yang R2-nya 78%)
RSS_Global <- sum(summary(rem_gls)$residuals^2) # RSS dari REM Two Ways
RSS_GWPR <- hasil.GWPR.tricube.fix$GW.diagnostic$RSS.gw # RSS dari GWPR Bisquare Fixed
# Ambil Degrees of Freedom (DF)
# Ganti obj_global dan obj_gwpr dengan nama objek model asli Anda
DF_Global <- rem_gls$df.residual
DF_GWPR <- hasil.GWPR.tricube.fix$GW.diagnostic$edf
# 2. Perhitungan Statistik F
# ------------------------------------------------------------------------------
# Hitung Selisih (Improvement)
SS_Improvement <- RSS_Global - RSS_GWPR # Seberapa banyak error berkurang
DF_Improvement <- DF_Global - DF_GWPR # Selisih derajat bebas
# Hitung Mean Square (Rata-rata Kuadrat)
MS_Improvement <- SS_Improvement / DF_Improvement
MS_GWPR <- RSS_GWPR / DF_GWPR
# Hitung F-Statistik
F_Hitung <- MS_Improvement / MS_GWPR
# Hitung P-Value
P_Value <- pf(F_Hitung, df1 = DF_Improvement, df2 = DF_GWPR, lower.tail = FALSE)
# 3. Menyusun Tabel Sesuai Format yang Anda Inginkan (F di Bawah)
# ------------------------------------------------------------------------------
anova_table <- data.frame(
Sumber_Variasi = c("Global Residuals (REM)", "GWPR Improvement", "GWPR Residuals"),
DF = c(DF_Global, DF_Improvement, DF_GWPR),
SS = c(RSS_Global, SS_Improvement, RSS_GWPR),
MS = c(NA, MS_Improvement, MS_GWPR),
F_Stat = c(NA, NA, F_Hitung),
P_Value = c(NA, NA, P_Value)
)
# Merapikan Angka (Pembulatan)
anova_table$SS <- round(anova_table$SS, 2)
anova_table$MS <- round(anova_table$MS, 2)
anova_table$F_Stat <- round(anova_table$F_Stat, 4)
anova_table$P_Value <- format(anova_table$P_Value, scientific = TRUE, digits = 4)
# Ubah NA menjadi kosong agar rapi saat diprint
anova_table[is.na(anova_table)] <- "NA"
# Tampilkan Hasil
print(anova_table)## Sumber_Variasi DF SS MS F_Stat P_Value
## 1 Global Residuals (REM) 144.00000 31795.53 NA NA NA
## 2 GWPR Improvement 36.07932 18451.80 511.42 NA NA
## 3 GWPR Residuals 107.92068 13343.73 123.64 4.1363 5.858e-0913. Visualisasi Peta Signifikansi
Membuat peta tematik yang menunjukkan kombinasi variabel yang signifikan berpengaruh terhadap kejadian banjir di setiap provinsi.
# VISUALISASI PETA: KOMBINASI VARIABEL SIGNIFIKAN ==============================
# 1. Persiapan Data Signifikansi
# Ambil hasil p-value yang sudah di-adjust. Pastikan urutan data sesuai.
# Karena data panel, kita ambil representasi unik provinsi (10 baris pertama)
# Asumsi: 10 baris pertama mewakili urutan provinsi di data 'banjir'
res_t_adj <- gwr.t.adjust(hasil.GWPR.tricube.fix)
p_vals_matrix <- res_t_adj$results$p[1:10, ]
# Buat Data Frame Signifikansi
df_sig_map <- data.frame(
PROV = banjir$PROV[1:10], # Mengambil nama provinsi
p_vals_matrix
)
# Tentukan nama variabel X
# (sesuaikan dengan colnames dari p_vals_matrix selain Intercept)
# Biasanya: X1.1, X2.1 dst atau X1, X2. Kita bersihkan nama kolomnya.
vars_target <- c("X1", "X2", "X3", "X4", "X5")
# Cek kolom mana yang sesuai dengan vars_target di df_sig_map
# (Terkadang output gwr memberi suffix .1, .2 dst)
cols_idx <- grep(paste(vars_target, collapse="|"), colnames(df_sig_map))
# 2. Membuat Kolom Kombinasi "var_sig"
df_sig_map$var_sig <- apply(df_sig_map[, cols_idx], 1, function(row) {
# Cek variabel mana yang p-value <= 0.05
is_sig <- row <= 0.05
# Ambil nama variabel yang TRUE
# Kita mapping index row ke nama vars_target agar rapi
nama_sig <- vars_target[is_sig]
if (length(nama_sig) == 0) {
return("Tidak Signifikan")
} else {
return(paste(sort(nama_sig), collapse = ", "))
}
})
# Cek tabel rekapitulasi sementara
print(table(df_sig_map$var_sig))##
## Tidak Signifikan X1, X5 X2, X3, X5 X2, X4
## 2 1 1 1
## X3 X3, X4, X5
## 2 3# 3. Persiapan Data Spasial (Shapefile)
library(rnaturalearth)
library(sf)
# Ambil peta Indonesia level provinsi
indo_map <- ne_states(country = "indonesia", returnclass = "sf")
# Filter hanya Pulau Sumatera
# Daftar nama provinsi di rnaturalearth mungkin Title Case, data Anda UPPERCASE.
# Kita samakan dulu.
indo_map$name_upper <- toupper(indo_map$name)
list_sumatera <- c("ACEH", "SUMATERA UTARA", "SUMATERA BARAT", "RIAU",
"KEPULAUAN RIAU", "JAMBI", "BENGKULU", "SUMATERA SELATAN",
"KEPULAUAN BANGKA BELITUNG", "LAMPUNG")
sumatera_sf <- indo_map %>%
filter(name_upper %in% list_sumatera)
# Gabungkan (Join) Data Signifikansi ke Peta
# Pastikan nama kolom kunci sama (PROV vs name_upper)
map_final <- sumatera_sf %>%
left_join(df_sig_map, by = c("name_upper" = "PROV"))
urutan_custom <- c("X2, X3, X5",
"X3, X4, X5",
"X1, X5",
"X2, X4",
"X3"
)
# Ubah var_sig menjadi factor agar legend rapi
map_final$var_sig <- factor(map_final$var_sig, levels = urutan_custom)
# 4. Visualisasi (Gaya Kustom)
# Definisikan palet warna dinamis
# kita pakai palette RColorBrewer
library(RColorBrewer)
jml_kategori <- length(levels(map_final$var_sig))
# Jika kategori sedikit pakai manual, jika banyak pakai palette otomatis
# 3. Definisikan Warna Berdasarkan Logika Pencampuran
# Format: c("KATEGORI" = "KODE_HEX")
# Definisikan warna pilihan
# Menambahkan warna Merah Tua (#800026) di bagian akhir
my_colors_logic <- c(
"#08519C",
"#3182BD",
"#9ECAE1",
"#31A354",
"#756BB1"
)
# Plotting
ggplot(data = map_final) +
geom_sf(aes(fill = var_sig), color = "black", size = 0.2) +
scale_fill_manual(
values = my_colors_logic,
name = "Peubah Signifikan",
na.value = "grey80", # Warna jika ada provinsi yang datanya NA
guide = guide_legend(
title.position = "top",
title.hjust = 0.5,
label.position = "right"
)
) +
labs(
title = "Peta Sebaran Signifikansi Parsial GWPR",
subtitle = "Wilayah Pulau Sumatera (2010-2024)",
caption = "Sumber: Hasil Olah Data (2026)"
) +
theme_void() + # Tema kosong (tanpa axis/grid)
theme(
legend.position = "right",
legend.title = element_text(face = "bold", size = 10),
legend.text = element_text(size = 9),
legend.key.width = unit(1, "cm"),
legend.key.height = unit(1, "cm"),
plot.title = element_text(hjust = 0.5, face = "bold",
size = 16, margin = margin(b=10)),
plot.subtitle = element_text(hjust = 0.5, size = 12, margin = margin(b=20))
)