1. Persiapan Pustaka (Library)

Tahap awal ini bertujuan untuk memuat seluruh paket (library) R yang diperlukan untuk manajemen data, analisis regresi data panel, pemodelan spasial (GWPR), hingga visualisasi pemetaan wilayah Sumatera.

2. Memuat dan Menyiapkan Data

Tahap ini dilakukan untuk menentukan direktori kerja, mengimpor dataset utama dari file Excel, serta melakukan penyesuaian nama kolom dan tipe data agar konsisten untuk analisis selanjutnya.

# 2. Data ======================================================================
setwd("D:/Kuliah/Lomba/PKM-AI 2026")
banjir <- read_excel("data_banjir_pulausumatera.xlsx", sheet = "2010")

# Penamaan Kolom
colnames(banjir) <- c("PROV", "TAHUN", "Y", "X1", "X2", "X3", "X4", "X5")
banjir$PROV <- as.factor(toupper(banjir$PROV))
banjir$TAHUN <- as.factor(banjir$TAHUN)

3. Pra-pemrosesan Data

Langkah ini mencakup prosedur pemeriksaan kualitas data melalui identifikasi nilai yang hilang (missing values), deteksi duplikasi data, serta visualisasi menggunakan boxplot untuk mengidentifikasi keberadaan outlier pada variabel dependen dan independen.

# 3. Pra-pemrosesan data =======================================================

# Cek Missing Value
print("Jumlah Missing Value:")

## [1] "Jumlah Missing Value:"

print(colSums(is.na(banjir)))

##  PROV TAHUN     Y    X1    X2    X3    X4    X5 
##     0     0     0     0     0     0     0     0

# Cek Redundansi Data
jumlah_duplikasi <- sum(duplicated(banjir))
print(paste("Jumlah data duplikat:", jumlah_duplikasi))

## [1] "Jumlah data duplikat: 0"

# Cek Outlier
df_long <- pivot_longer(banjir, cols = c("Y", "X1", "X2", "X3", "X4", "X5"), 
                        names_to = "Variabel", values_to = "Nilai")

df_long$Variabel <- factor(df_long$Variabel, 
                           levels = c("Y", "X1", "X2", "X3", "X4", "X5"))

p_clean <- ggplot(df_long, aes(x = Variabel, y = Nilai)) +
  geom_boxplot(fill = "grey95", 
               color = "black",         
               lwd = 0.6,                
               outlier.colour = "red",  
               outlier.shape = 19,        
               outlier.size = 1.5,
               outlier.alpha = 0.6) +       
scale_y_continuous(
  labels = scales::label_number(big.mark = ".", decimal.mark = ",")
) + 
facet_wrap(~Variabel, scales = "free", ncol = 6) + 
  labs(title = "Distribusi Data dan Deteksi Outlier",
       y = "", 
       x = "") + 
  theme_bw() +
  theme(
    text = element_text(family = "serif"), 
    plot.title = element_text(face = "bold", size = 14, hjust = 0),
    plot.subtitle = element_text(face = "italic", size = 11, hjust = 0, color = "grey30"),
    strip.background = element_rect(fill = "white", color = "black"), 
    strip.text = element_text(face = "bold", size = 11),
    axis.text = element_text(color = "black"),
    panel.grid.minor = element_blank(), 
    panel.grid.major.x = element_blank() 
  )
print(p_clean)

4. Analisis Statistik Deskriptif

Bagian ini menyajikan ringkasan statistik deskriptif untuk seluruh variabel penelitian, yang meliputi nilai minimum, maksimum, rata-rata (mean), dan standar deviasi guna memahami karakteristik distribusi data secara umum.

# 4. Analisis statistik deskriptif =============================================
summary_stats <- banjir %>%
  dplyr::select("Y", "X1", "X2", "X3", "X4","X5") %>%
  summarise_all(list(
    Min = ~min(.),
    Max = ~max(.),
    Mean = ~mean(.),
    SD = ~sd(.)
  )) %>%
  pivot_longer(everything(),
               names_to = c("Variable", "Statistic"),
               names_sep = "_") %>%
  pivot_wider(names_from = Statistic, values_from = value)
kable(summary_stats, caption = "Tabel 2. Statistik deskriptif variabel penelitian")

5. Seleksi Peubah dan Uji Multikolinearitas

Langkah ini dilakukan untuk memastikan tidak adanya hubungan linear yang kuat antar variabel independen melalui pengujian Variance Inflation Factor (VIF). Hal ini krusial agar estimasi parameter dalam model GWPR tetap objektif dan terhindar dari bias standar error.

# 5. Seleksi Peubah Multikolinearitas ==========================================
check_model <- lm(Y ~ X1 + X2 + X3 + X4 + X5, banjir)
multikol <- data.frame(t(vif(check_model)))
kable(multikol, caption = "Tabel 3. Hasil uji multikolinearitas menggunakan VIF ")

6. Estimasi Common Effect Model (CEM)

Tahap ini melakukan pemodelan awal menggunakan pendekatan Common Effect Model (CEM) atau pooling regression. Model ini mengasumsikan bahwa perilaku data antar individu (provinsi) dan antar waktu adalah sama, sehingga estimasi dilakukan dengan menggabungkan seluruh data tanpa memperhatikan efek spesifik individu maupun waktu.

# 6. 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

7. Estimasi Fixed Effect Model (FEM)

Pada bagian ini, dilakukan estimasi menggunakan pendekatan Fixed Effect Model (FEM) dengan mempertimbangkan efek spesifik yang tidak tertangkap oleh variabel independen. Pemodelan mencakup tiga skenario: Individual Effect (efek antar provinsi), Time Effect (efek antar tahun), dan Two-Ways Effect (gabungan keduanya). Dilakukan juga uji Breusch-Pagan (BP Test) untuk mendeteksi keberadaan efek galat yang tidak teramati (unobserved effects) pada model.

# 7. 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-11

summary(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.02875

summary(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.025186

summary(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.9905839

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 effects

plmtest(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 effects

plmtest(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

8. Estimasi Random Effect Model (REM)

Tahap ini melakukan estimasi menggunakan pendekatan Random Effect Model (REM) dengan metode Nerlove. Berbeda dengan FEM, model ini mengasumsikan bahwa perbedaan antar individu bersifat acak dan tidak berkorelasi dengan variabel independen. Selanjutnya, dilakukan uji Lagrange Multiplier (Breusch-Pagan) untuk membandingkan kebaikan model REM terhadap CEM berdasarkan efek individu, waktu, maupun keduanya.

# 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

9. Uji Pemilihan Model dan Tabel Rekapitulasi

Bagian ini merupakan tahap krusial untuk menentukan model data panel yang paling tepat melalui tiga pengujian statistik: Uji Chow (CEM vs FEM), Uji Hausman (FEM vs REM), dan Uji Lagrange Multiplier (CEM vs REM). Seluruh hasil pengujian dirangkum dalam tabel rekapitulasi otomatis untuk mempermudah pengambilan keputusan model terbaik yang akan digunakan dalam analisis lebih lanjut.

# 9. 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, "FEM lebih baik dari CEM", "CEM lebih baik dari FEM"))
  } else if (test_type == "Hausman") {
    return(ifelse(p_val < alpha, "FEM lebih baik dari REM", "REM lebih baik dari FEM"))
  } else if (test_type == "LM") {
    return(ifelse(p_val < alpha, "REM lebih baik dari CEM", "CEM lebih baik dari REM"))
  }
}

# --- MEMBUAT DATA FRAME HASIL ---
hasil_uji <- data.frame(
  Uji = c("Chow", "LM", "Hausman"),
  
  Statistik = c(
    as.numeric(chow_test$statistic), 
    as.numeric(lm_test$statistic),
    as.numeric(hausman_test$statistic)
  ),
  
  P_Value = c(
    as.numeric(chow_test$p.value), 
    as.numeric(lm_test$p.value),
    as.numeric(hausman_test$p.value)
  ),
  
  Kesimpulan = c(
    get_conclusion(as.numeric(chow_test$p.value), "Chow"),
    get_conclusion(as.numeric(lm_test$p.value), "LM"),
    get_conclusion(as.numeric(hausman_test$p.value), "Hausman")
  )
)

# --- FORMATTING TAMPILAN AGAR RAPI (OPSIONAL) ---
# Membulatkan angka agar enak dibaca
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)
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, caption = "Tabel 4. Hasil Pemilihan Model Regresi Data Panel")

10. Uji Diagnostik Sisaan (Asumsi Klasik dan Dependensi Spasial)

Tahap ini bertujuan untuk memvalidasi kelayakan model Random Effect Model (REM) melalui serangkaian uji asumsi klasik yang meliputi Normalitas, Heteroskedastisitas, dan Autokorelasi. Selain itu, dilakukan pula uji Pesaran CD (Cross-sectional Dependence) untuk mendeteksi adanya ketergantungan spasial antar wilayah. Hasil dari tahap ini akan memberikan gambaran apakah model panel global sudah memadai atau diperlukan pendekatan lokal seperti GWPR.

# 10. 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)
uji_auto_serial <- pbgtest(rem_gls)
# D. UJI BARU: Uji Dependensi Spasial / Cross-sectional Dependence (Pesaran CD)
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)", 
                 "Homoskedastis (BP)", 
                 "Autokorelasi (BG)", 
                 "Dependensi Spasial (CD)"), # Baris Baru
  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)),
  Statistik = c(round(uji_norm$statistic, 4), 
                round(uji_hetero$statistic, 4), 
                round(uji_auto_serial$statistic, 4), 
                round(uji_spasial$statistic, 4)),
  Keputusan = c(ifelse(uji_norm$p.value > 0.05,
                       "Terpenuhi", "Tidak Terpenuhi"),
                ifelse(uji_hetero$p.value > 0.05,
                       "Terpenuhi", "Tidak Terpenuhi"),
                ifelse(uji_auto_serial$p.value < 0.05,
                       "Terpenuhi", "Tidak Terpenuhi"),
                ifelse(uji_spasial$p.value < 0.05,
                       "Terpenuhi", "Tidak Terpenuhi"))
)
# --- 4. Cetak Output dengan Kable ---
kable(tabel_estimasi, caption = "Tabel 5.1. Koefisien Estimasi REM dan Uji Multikolinearitas")

kable(tabel_diagnostik, caption = "Tabel 5.2. Hasil diagnostik asumsi klasik dan efek spasial")

11. Geographically Weighted Panel Regression (GWPR): Integrasi Data Spasial

Pada tahap ini, dilakukan persiapan data untuk pemodelan GWPR dengan mengintegrasikan koordinat geografis (Latitude dan Longitude) ke dalam dataset panel. Proses ini melibatkan penggabungan data tabular dengan referensi spasial wilayah provinsi di Pulau Sumatera serta transformasi objek data menjadi format spasial (SpatialPointsDataFrame) agar dapat diproses dalam estimasi regresi berbobot geografis.

# 11. 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"

12. Pemilihan Fungsi Pembobot Spasial Terbaik

Tahap ini bertujuan untuk menentukan fungsi kernel dan jenis bandwidth yang paling optimal untuk model GWPR. Pengujian dilakukan dengan membandingkan dua kelompok fungsi pembobot, yaitu Adaptive Kernel (jarak fleksibel berdasarkan kepadatan data) dan Fixed Kernel (jarak tetap), menggunakan empat kriteria fungsi: Bisquare, Gaussian, Exponential, dan Tricube. Penentuan model terbaik didasarkan pada nilai AICc terkecil, R-Squared Adjusted tertinggi, serta Residual Sum of Squares (RSS) minimum guna menjamin akurasi estimasi lokal.

# 12. Pemilihan Fungsi Pembobot Terbaik ========================================
# --- 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.15

hasil.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.64

hasil.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.57

hasil.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.9

hasil.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.33

hasil.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.98

hasil.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.69

hasil.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.05

hasil.GWPR.tricube.fix <- gwr.basic(Y ~ X1 + X2 + X3 + X4 + X5,
                                    data = data.sp.GWPR,
                                    bw = bwd.GWPR.tricube.fix,
                                    kernel = "tricube", adaptive=F)


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)   

# Gabungkan menjadi tabel
tabel <- cbind(RSS, R2adj, AICc)

# Beri nama baris untuk 8 model
rownames(tabel) <- c("Adaptive Bisquare", 
                     "Adaptive Gaussian",
                     "Adaptive Exponential", 
                     "Adaptive Tricube", 
                     "Fixed Bisquare",
                     "Fixed Gaussian", 
                     "Fixed Exponential",
                     "Fixed Tricube")

# Tampilkan hasil
kable(tabel, caption = "Tabel 6. Evaluasi Kriteria Statistik untuk Pemilihan Fungsi Pembobot Terbaik")

13. Uji ANOVA: Perbandingan Model Global (REM) vs Model Lokal (GWPR)

Tahap ini bertujuan untuk memvalidasi secara statistik apakah penggunaan model GWPR memberikan peningkatan kualitas model yang signifikan dibandingkan dengan model global (Random Effect Model). Melalui perbandingan nilai Residual Sum of Squares (RSS) dan Adjusted R-Square, uji ANOVA ini menghitung nilai F-Statistik untuk membuktikan apakah keberadaan variasi spasial dalam data determinan banjir di Pulau Sumatera memang memerlukan pendekatan model lokal.

# 13. Uji ANOVA (UJI F) PERBANDINGAN MODEL GLOBAL VS GWPR
R2adj_Global <- summary(rem_gls)$r.squared[2]
R2adj_GWPR   <- hasil.GWPR.tricube.fix$GW.diagnostic$gwR2.adj

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

DF_Global  <- rem_gls$df.residual               
DF_GWPR    <- hasil.GWPR.tricube.fix$GW.diagnostic$edf 

R2adj_Improvement <- R2adj_Global - R2adj_GWPR       # Seberapa banyak error berkurang
SS_Improvement    <- RSS_Global - RSS_GWPR       # Seberapa banyak error berkurang
DF_Improvement    <- DF_Global - DF_GWPR         # Selisih derajat bebas

MS_Improvement <- SS_Improvement / DF_Improvement
MS_GWPR        <- RSS_GWPR / DF_GWPR

F_Hitung <- MS_Improvement / MS_GWPR

P_Value  <- pf(F_Hitung, df1 = DF_Improvement, df2 = DF_GWPR, lower.tail = FALSE)

anova_table <- data.frame(
  Sumber_Variasi = c("Global Residuals (REM)", "GWPR Improvement", "GWPR Residuals"),
  
  DF = c(DF_Global, DF_Improvement, DF_GWPR),
  
  R2adj = c(R2adj_Global, R2adj_Improvement, R2adj_GWPR),
  
  F_Stat = c(NA, NA, F_Hitung),
  
  P_Value = c(NA, NA, P_Value)
)

# Merapikan Angka (Pembulatan)
anova_table$R2adj <- abs(round(anova_table$R2adj, 4))
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# Taabs()mpilkan Hasil
kable(anova_table, caption = "Tabel 7. Hasil Uji ANOVA perbandingan kinerja Model Global dan GWPR")

14. Uji t dan Tipologi Signifikansi Parameter Lokal Model GWPR

Tahap ini dilakukan untuk menguji signifikansi parameter secara lokal pada setiap lokasi pengamatan di Pulau Sumatera. Melalui Uji t yang telah disesuaikan (t-adjust), kita dapat mengidentifikasi variabel independen mana saja yang berpengaruh signifikan terhadap kejadian banjir di tiap provinsi. Hasil pengujian ini kemudian digunakan untuk menyusun Tipologi Wilayah, yaitu pengelompokan provinsi berdasarkan kombinasi peubah signifikan yang serupa, guna memberikan rekomendasi kebijakan mitigasi bencana yang lebih spesifik dan tepat sasaran di tingkat lokal.

# 14. Uji t (Distribusi Spasial Signifikansi Parameter Lokal Model GWPR) =======
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,]

df_gabung <- p_value %>%
  mutate(PROV = banjir$PROV[1:10]) 

tabel_final <- df_gabung %>%
  # Ubah data dari format lebar (wide) ke panjang (long)
  pivot_longer(
    cols = -PROV,             # Pilih semua kolom kecuali PROV
    names_to = "Variabel",    # Nama kolom baru untuk X1, X2, dst
    values_to = "Status"      # Isi nilainya (Signifikan/Tidak)
  ) %>%
  # Filter hanya ambil yang statusnya "Signifikan"
  filter(Status == "Signifikan") %>%
  # Kelompokkan berdasarkan nama Variabel
  group_by(Variabel) %>%
  # Lakukan rekapitulasi
  summarise(
    Daftar_Provinsi = paste(PROV, collapse = ", "), # Menggabungkan nama provinsi dengan koma
    Jumlah = n()                                    # Menghitung jumlah provinsi
  ) %>%
  # (Opsional) Bersihkan nama variabel jika ada suffix "_p" agar lebih rapi
  mutate(Variabel = gsub("_p", "", Variabel))

kable(tabel_final, caption = "Tabel 8. Distribusi Spasial Signifikansi Parameter Lokal Model GWPR ")

# 1. Persiapan Data Signifikansi 
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
)
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

tabel_tipologi <- df_sig_map %>%
  group_by(var_sig) %>%
  summarise(
    daftar_prov = paste(PROV, collapse = ", "),
    .groups = 'drop'
  ) %>%
  mutate(Kelompok = row_number()) %>%
  select(
    Kelompok, 
    `Peubah Signifikan` = var_sig, 
    `Daftar Provinsi` = daftar_prov
  )

kable(tabel_tipologi, caption = "Tabel 9. Pengelompokkan wilayah berdasarkan peubah signifikan pada Model GWPR")

15. Distribusi Spasial Determinan Banjir di Pulau Sumatera Periode 2010–2024

Bagian ini menyajikan visualisasi geografis dari hasil pemodelan GWPR menggunakan peta tematik. Visualisasi dilakukan dengan memetakan pengelompokan provinsi berdasarkan variabel-variabel yang berpengaruh signifikan di masing-masing wilayah. Melalui representasi spasial ini, perbedaan karakteristik pendorong banjir antar wilayah di Pulau Sumatera dapat diidentifikasi secara visual, sehingga mempermudah interpretasi terhadap heterogenitas spasial yang terjadi selama periode pengamatan.

# 15. Distribusi Spasial Determinan Banjir di Pulau Sumatera Periode 2010–2024


indo_map <- ne_states(country = "indonesia", returnclass = "sf")

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 
jml_kategori <- length(levels(map_final$var_sig))

my_colors_logic <- c(
  "#000080",  # 1. Biru Tua (Navy)
  "#660099",  # 2. Ungu (Purple)
  "#FF1493",  # 3. Merah Muda (Deep Pink/Magenta)
  "#FF8C00",  # 4. Oranye (Dark Orange)
  "#FFFF00"   # 5. Kuning (Yellow)
)
# 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: TIM PKM-AI Statistika 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))
  )