Disclaimer: Modul pembelajaran ini disusun dengan menyintesis berbagai sumber literatur ilmiah yang tercantum pada bagian referensi. Dalam proses pengerjaannya, teknologi kecerdasan buatan (Artificial Intelligence/AI) digunakan sebagai alat bantu untuk meringkas, menyatukan informasi dari beberapa sumber sekaligus, serta menyusun draf sintaks pemrograman. Kendati demikian, seluruh isi, logika pemodelan, dan kode pemrograman di dalam modul ini telah melalui proses pengecekan, pengujian, serta validasi secara mandiri oleh penyusun sebelum diterbitkan.

Pengantar: Kapan dan Mengapa GWR?

Latar Belakang

Dalam analisis spasial, kita sering berasumsi bahwa hubungan antarvariabel bersifat stasioner — yaitu konstan di seluruh wilayah studi. Regresi Kuadrat Terkecil Biasa (Ordinary Least Squares/OLS) mengikuti asumsi ini dengan menghasilkan satu set koefisien global yang mewakili rata-rata hubungan di seluruh area.

Namun, kenyataan geografis menunjukkan hal berbeda. Tobler (1970) menyatakan dalam Hukum Pertama Geografi: “Everything is related to everything else, but near things are more related than distant things.” Prinsip ini menyiratkan bahwa hubungan antarvariabel dapat berbeda-beda bergantung pada lokasi.

Sebagai contoh:

  • Pengaruh pendapatan terhadap harga rumah mungkin lebih kuat di pusat kota daripada di pinggiran.

  • Pengaruh kemiskinan terhadap angka kriminalitas mungkin berbeda antara wilayah perkotaan dan perdesaan.

  • Pengaruh elevasi terhadap curah hujan bervariasi di setiap wilayah topografi.

Fenomena ini disebut non-stasioneritas spasial atau heterogenitas spasial dalam hubungan antar-variabel.

Apa itu GWR?

Geographically Weighted Regression (GWR) adalah teknik regresi yang memungkinkan koefisien regresi bervariasi secara spasial. GWR dikembangkan oleh Brunsdon, Fotheringham, dan Charlton (1996) sebagai ekstensi lokal dari model regresi global.

Alih-alih menghasilkan satu set koefisien global, GWR menghasilkan satu set koefisien unik untuk setiap lokasi dalam dataset. Estimasi dilakukan dengan memberikan bobot lebih besar pada pengamatan yang lebih dekat secara geografis dengan lokasi yang sedang diestimasi, mirip dengan prinsip kernel density estimation.

Kapan Menggunakan GWR?

GWR tepat digunakan ketika:

Kondisi Indikator
Dugaan non-stasioneritas spasial Ada alasan teoritis mengapa hubungan antar-variabel bervariasi antar-lokasi
Residual OLS berautokorelasi spasial Moran’s I dari residual OLS signifikan
Heterogenitas spasial terdeteksi Uji Breusch-Pagan atau Koenker-Bassett signifikan
Tujuan eksplorasi spasial Ingin memetakan di mana suatu variabel memiliki pengaruh positif vs. negatif
Data bersifat agregat spasial Data pada level wilayah (kecamatan, kabupaten, dll.)

GWR kurang tepat digunakan ketika: - Data tidak memiliki komponen spasial yang bermakna - Jumlah pengamatan sangat sedikit (n < 50) - Tujuan utama adalah prediksi bukan eksplorasi (GWR cenderung overfitting) - Multikolinieritas lokal sangat tinggi

Perbedaan OLS vs. GWR

Aspek OLS (Regresi Global) GWR (Regresi Lokal)
Koefisien Satu nilai untuk seluruh wilayah Berbeda di setiap lokasi
Asumsi Stasioneritas spasial Non-stasioneritas diperbolehkan
Output Tabel koefisien Peta koefisien lokal
Interpretasi Rata-rata pengaruh global Pengaruh yang bervariasi secara spasial
Kompleksitas Sederhana Lebih kompleks
Overfitting Rendah Perlu kehati-hatian

Konsep Pemodelan GWR

Model OLS (Regresi Global)

Model regresi OLS standar:

\[y_i = \beta_0 + \sum_{k=1}^{p} \beta_k x_{ik} + \varepsilon_i\]

di mana \(\beta_0, \beta_1, \ldots, \beta_p\) adalah koefisien global (sama untuk semua \(i\)), dan \(\varepsilon_i \sim N(0, \sigma^2)\).

Estimasi OLS: \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\)

2.2 Model GWR

Model GWR memperluas OLS dengan membiarkan koefisien bervariasi secara spasial:

\[y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i)\, x_{ik} + \varepsilon_i\]

di mana \((u_i, v_i)\) adalah koordinat geografis lokasi ke-\(i\), dan \(\beta_k(u_i, v_i)\) adalah koefisien lokal yang bergantung pada lokasi.

Estimasi Parameter GWR

Koefisien GWR diestimasi dengan Weighted Least Squares (WLS) secara terpisah untuk setiap lokasi:

\[\hat{\boldsymbol{\beta}}(u_i, v_i) = \left[\mathbf{X}^T \mathbf{W}(u_i, v_i)\, \mathbf{X}\right]^{-1} \mathbf{X}^T \mathbf{W}(u_i, v_i)\, \mathbf{y}\]

di mana \(\mathbf{W}(u_i, v_i)\) adalah matriks bobot diagonal \(n \times n\):

\[\mathbf{W}(u_i, v_i) = \mathrm{diag}\left[w_{i1},\, w_{i2},\, \ldots,\, w_{in}\right]\]

Bobot \(w_{ij}\) mencerminkan kedekatan spasial antara lokasi ke-\(i\) (lokasi yang sedang diestimasi) dan lokasi ke-\(j\) (pengamatan yang digunakan sebagai data). Pengamatan yang lebih dekat mendapat bobot lebih besar.

Fungsi Kernel

Fungsi kernel menentukan bagaimana bobot menurun seiring bertambahnya jarak \(d_{ij}\) dari titik estimasi.

a. Gaussian (kontinu, batas tak terhingga): \[w_{ij} = \exp\!\left(-\frac{d_{ij}^2}{2h^2}\right)\]

b. Exponential (kontinu, batas tak terhingga): \[w_{ij} = \exp\!\left(-\frac{d_{ij}}{h}\right)\]

c. Bisquare (diskontinu, batas terhingga — compact support): \[w_{ij} = \begin{cases} \left[1 - \left(\dfrac{d_{ij}}{h}\right)^2\right]^2 & \text{jika } d_{ij} \leq h \\ 0 & \text{jika } d_{ij} > h \end{cases}\]

d. Tricube (diskontinu, compact support): \[w_{ij} = \begin{cases} \left[1 - \left(\dfrac{d_{ij}}{h}\right)^3\right]^3 & \text{jika } d_{ij} \leq h \\ 0 & \text{jika } d_{ij} > h \end{cases}\]

e. Boxcar (semua pengamatan dalam radius mendapat bobot sama): \[w_{ij} = \begin{cases} 1 & \text{jika } d_{ij} \leq h \\ 0 & \text{jika } d_{ij} > h \end{cases}\]

Catatan: Kernel bisquare dan tricube memiliki compact support — pengamatan di luar bandwidth mendapat bobot nol. Ini membuat komputasi lebih efisien. Kernel Gaussian dan exponential secara teoritis menggunakan semua pengamatan, tetapi pengamatan jauh hampir tidak berpengaruh.

enis Bandwidth

Ada dua jenis bandwidth \(h\):

a. Fixed bandwidth (bandwidth tetap):
Jari-jari pencarian konstan dalam satuan jarak (km, meter). Cocok ketika data terdistribusi merata secara spasial. Kelemahan: di daerah data padat, banyak pengamatan digunakan; di daerah jarang, sedikit pengamatan tersedia.

b. Adaptive bandwidth (bandwidth adaptif):
Jumlah tetangga (k nearest neighbors) konstan di setiap lokasi. Bandwidth (jarak) menyesuaikan kepadatan data. Di daerah padat, bandwidth sempit; di daerah jarang, bandwidth lebar.

Ukuran Kebaikan Model GWR

AICc (Corrected Akaike Information Criterion): \[\text{AICc} = 2n\ln(\hat{\sigma}) + n\ln(2\pi) + n\left(\frac{n + \mathrm{tr}(\mathbf{H})}{n - 2 - \mathrm{tr}(\mathbf{H})}\right)\]

R² Lokal: Mengukur kecocokan model di setiap lokasi.
R² Global GWR: Rata-rata tertimbang R² lokal.

Persiapan: Paket dan Data

Instalasi dan Pemuatan Paket

# Instalasi paket (jalankan sekali saja)
install.packages(c(
  "GWmodel",    # GWR utama
  "sf",         # simple features (data spasial modern)
  "spdep",      # uji dependensi spasial
  "tmap",       # visualisasi peta tematik
  "ggplot2",    # visualisasi umum
  "tidyverse",  # manipulasi data
  "car",        # uji regresi (VIF, dll.)
  "lmtest",     # uji diagnostik regresi
  "RColorBrewer"# palet warna
))
# Memuat paket
library(GWmodel)
library(sf)
library(spdep)
library(tmap)
library(ggplot2)
library(tidyverse)
library(car)
library(lmtest)
library(RColorBrewer)

# Pengaturan tampilan peta
tmap_mode("plot")

Dataset: Georgia County-Level Data

Dataset yang digunakan adalah Georgia, yang tersedia langsung dalam paket GWmodel. Dataset ini berisi data tingkat kabupaten (county) untuk negara bagian Georgia, Amerika Serikat, tahun 1990 (n = 159 kabupaten).

# Memuat data Georgia dari paket GWmodel
data(Georgia) #akan tersimpan sebagai Gedu.df

Georgia<-Gedu.df

# Melihat kelas objek
class(Georgia)
## [1] "data.frame"
# Melihat nama variabel
names(Georgia)
##  [1] "AreaKey"  "Latitude" "Longitud" "TotPop90" "PctRural" "PctBach" 
##  [7] "PctEld"   "PctFB"    "PctPov"   "PctBlack" "ID"       "X"       
## [13] "Y"

Keterangan variabel:

Variabel Deskripsi
AreaKey Nomor identifikasi untuk setiap county (kabupaten/wilayah administratif)
Latitude Garis lintang (latitude) titik pusat geografis county
Longitud Garis bujur (longitude) titik pusat geografis county
TotPop90 Jumlah penduduk county pada tahun 1990
PctRural Persentase penduduk county yang dikategorikan sebagai penduduk pedesaan (%)
PctBach Persentase penduduk county yang memiliki gelar sarjana (%)
PctEld Persentase penduduk county yang berusia 65 tahun atau lebih (%)
PctFB Persentase penduduk county yang lahir di luar Amerika Serikat (%)
PctPov Persentase penduduk county yang hidup di bawah garis kemiskinan (%)
PctBlack Persentase penduduk county yang merupakan kelompok ras kulit hitam (%)

Variabel respons: PctBach (Persentase gelar sarjana)
Variabel prediktor: PctFB, PctPov, PctBlack, PctEld

Pertanyaan penelitian: Apakah faktor-faktor sosial ekonomi (kemiskinan, komposisi etnis, imigran, proporsi lansia) memengaruhi persentase penduduk berpendidikan sarjana, dan apakah pengaruh ini bervariasi secara spasial?

Import dan Eksplorasi Data

Konversi ke Format sf

# Konversi SpatialPolygonsDataFrame ke sf
georgia_sf <- st_as_sf(Georgia,coords = c("X", "Y"))

# Lihat struktur data
str(georgia_sf, max.level = 2)
## Classes 'sf' and 'data.frame':   159 obs. of  12 variables:
##  $ AreaKey : int  13001 13003 13005 13007 13009 13011 13013 13015 13017 13019 ...
##  $ Latitude: num  31.8 31.3 31.6 31.3 33.1 ...
##  $ Longitud: num  -82.3 -82.9 -82.5 -84.5 -83.3 ...
##  $ TotPop90: int  15744 6213 9566 3615 39530 10308 29721 55911 16245 14153 ...
##  $ PctRural: num  75.6 100 61.7 100 42.7 100 64.6 75.2 47 66.2 ...
##  $ PctBach : num  8.2 6.4 6.6 9.4 13.3 6.4 9.2 9 7.6 7.5 ...
##  $ PctEld  : num  11.43 11.77 11.11 13.17 8.64 ...
##  $ PctFB   : num  0.64 1.58 0.27 0.11 1.43 0.34 0.92 0.82 0.33 1.19 ...
##  $ PctPov  : num  19.9 26 24.1 24.8 17.5 15.1 14.7 10.7 22 19.3 ...
##  $ PctBlack: num  20.8 26.9 15.4 51.7 42.4 ...
##  $ ID      : int  133 158 146 155 79 23 33 24 138 153 ...
##  $ geometry:sfc_POINT of length 159; first list element:  'XY' num  941397 3521764
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:11] "AreaKey" "Latitude" "Longitud" "TotPop90" ...
# Lihat beberapa baris pertama
head(georgia_sf[, c("PctBach", "PctFB", "PctPov","PctBlack", "PctEld")], 10)
## Simple feature collection with 10 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 699011.5 ymin: 3466377 xmax: 941396.6 ymax: 3807616
## CRS:           NA
##    PctBach PctFB PctPov PctBlack PctEld                 geometry
## 1      8.2  0.64   19.9    20.76  11.43 POINT (941396.6 3521764)
## 2      6.4  1.58   26.0    26.86  11.77   POINT (895553 3471916)
## 3      6.6  0.27   24.1    15.42  11.11 POINT (930946.4 3502787)
## 4      9.4  0.11   24.8    51.67  13.17 POINT (745398.6 3474765)
## 5     13.3  1.43   17.5    42.39   8.64 POINT (849431.3 3665553)
## 6      6.4  0.34   15.1     3.49  11.37 POINT (819317.3 3807616)
## 7      9.2  0.92   14.7    11.44  10.63 POINT (803747.1 3769623)
## 8      9.0  0.82   10.7     9.21   9.66 POINT (699011.5 3793408)
## 9      7.6  0.33   22.0    31.33  12.81 POINT (863020.8 3520432)
## 10     7.5  1.19   19.3    11.62  11.98 POINT (859915.8 3466377)

Statistik Deskriptif

# Statistik deskriptif variabel utama
vars_interest <- c("PctBach", "PctFB", "PctPov", "PctBlack", "PctEld")

summary(st_drop_geometry(georgia_sf)[, vars_interest])
##     PctBach          PctFB           PctPov         PctBlack    
##  Min.   : 4.20   Min.   :0.040   Min.   : 2.60   Min.   : 0.00  
##  1st Qu.: 7.60   1st Qu.:0.415   1st Qu.:14.05   1st Qu.:11.75  
##  Median : 9.40   Median :0.720   Median :18.60   Median :27.64  
##  Mean   :10.95   Mean   :1.131   Mean   :19.34   Mean   :27.39  
##  3rd Qu.:12.00   3rd Qu.:1.265   3rd Qu.:24.65   3rd Qu.:40.06  
##  Max.   :37.50   Max.   :6.740   Max.   :35.90   Max.   :79.64  
##      PctEld     
##  Min.   : 1.46  
##  1st Qu.: 9.81  
##  Median :12.07  
##  Mean   :11.74  
##  3rd Qu.:13.70  
##  Max.   :22.96
# Distribusi variabel dengan histogram
georgia_sf %>% 
  st_drop_geometry() %>% 
  select(all_of(vars_interest)) %>% 
  pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai") %>% 
  ggplot(aes(x = Nilai)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white") +
  facet_wrap(~Variabel, scales = "free") +
  labs(title = "Distribusi Variabel Utama",
       subtitle = "Data Georgia, 1990",
       x = "Nilai", y = "Frekuensi") +
  theme_bw()

Matriks Korelasi

# Matriks korelasi antar variabel
cor_matrix <- cor(st_drop_geometry(georgia_sf)[, vars_interest], 
                  use = "complete.obs")
round(cor_matrix, 3)
##          PctBach  PctFB PctPov PctBlack PctEld
## PctBach    1.000  0.672 -0.402   -0.109 -0.458
## PctFB      0.672  1.000 -0.329   -0.112 -0.483
## PctPov    -0.402 -0.329  1.000    0.736  0.568
## PctBlack  -0.109 -0.112  0.736    1.000  0.297
## PctEld    -0.458 -0.483  0.568    0.297  1.000
# Visualisasi matriks korelasi
# install.packages("corrplot")  # jika belum terinstal
library(corrplot)
corrplot(cor_matrix, method = "color", type = "upper",
         tl.cex = 0.8, addCoef.col = "black", number.cex = 0.7,
         title = "Matriks Korelasi Variabel Utama",
         mar = c(0, 0, 1, 0))

Visualisasi Spasial Variabel Respons

data("GeorgiaCounties") #akan tersimpan sebagai Gedu.counties

georgia_poly <- st_as_sf(Gedu.counties)

# Lihat nama variabel pada kedua data
names(Georgia)
##  [1] "AreaKey"  "Latitude" "Longitud" "TotPop90" "PctRural" "PctBach" 
##  [7] "PctEld"   "PctFB"    "PctPov"   "PctBlack" "ID"       "X"       
## [13] "Y"
names(georgia_poly)
## [1] "AREA"      "PERIMETER" "G_UTM_"    "G_UTM_ID"  "AREANAME"  "AREAKEY"  
## [7] "X_COORD"   "Y_COORD"   "geometry"
class(georgia_poly$AREAKEY)
## [1] "factor"
class(Georgia$AreaKey)
## [1] "integer"
georgia_poly <- georgia_poly %>%
  mutate(AREAKEY = as.integer(as.character(AREAKEY)))

# Gabungkan atribut Georgia ke poligon
# (umumnya menggunakan AreaKey)
georgia_map <- georgia_poly %>%
  left_join(Georgia, by = c("AREAKEY" = "AreaKey"))

# Cek hasil join
summary(georgia_map$PctBach)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.20    7.60    9.40   10.95   12.00   37.50
# Peta kombinasi: Poligon latar belakang abu-abu + tumpukan titik besar (tmap v4)
map_kombinasi <- tm_shape(georgia_map) +
  tm_polygons(fill = "grey92", col = "grey70", lwd = 0.4) + # Background peta abu-abu
  
  tm_shape(georgia_sf) +
  tm_dots(
    fill = "PctBach",
    size = 0.6, # Mengunci ukuran titik menjadi lebih besar dan terlihat jelas
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "YlOrRd"),
    fill.legend = tm_legend(title = "PctBach (%)", position = tm_pos_out("right"))
  ) +
  tm_title("Georgia: Persentase Penduduk dengan Gelar Sarjana (1990)", size = 0.85) +
  tm_compass(position = tm_pos_in("left", "bottom")) +
  tm_scalebar(position = tm_pos_in("right", "bottom"))

map_kombinasi

# Buat peta
tm_shape(georgia_map) +
  tm_polygons(fill = "PctBach",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5,
      values = "YlOrRd"),
    fill.legend = tm_legend(title = "PctBach (%)"),
    col = "grey40",
    lwd = 0.4) +
  tm_title("Georgia: Persentase Penduduk dengan Gelar Sarjana (1990)") +
  tm_compass(position = c("left", "bottom")) +
  tm_scalebar(position = c("right", "bottom")
  )

# Peta distribusi variabel prediktor
st_crs(georgia_map) <- 26917

tm_shape(georgia_map) +
  tm_polygons(fill = "PctFB",
              fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.blues"),
              fill.legend = tm_legend(title = "Persentase (%)"),
              col = "grey60", lwd = 0.2) +
  tm_title("Georgia: Persentase Penduduk Lahir di Luar Amerika Serikat (1990)") +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("left", "bottom")) +
  tm_layout(legend.outside = TRUE, frame = FALSE)
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

tm_shape(georgia_map) +
  tm_polygons(fill = "PctPov",
              fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.reds"),
              fill.legend = tm_legend(title = "Persentase (%)"),
              col = "grey60", lwd = 0.2) +
  tm_title("Georgia: Persentase Penduduk di Bawah Garis Kemiskinan (1990)") +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("left", "bottom")) +
  tm_layout(legend.outside = TRUE, frame = FALSE)
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

tm_shape(georgia_map) +
  tm_polygons(fill = "PctBlack",
              fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.purples"),
              fill.legend = tm_legend(title = "Persentase (%)"),
              col = "grey60", lwd = 0.2) +
  tm_title("Georgia: Persentase Penduduk Black (1990)") +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("left", "bottom")) +
  tm_layout(legend.outside = TRUE, frame = FALSE)
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

tm_shape(georgia_map) +
  tm_polygons(fill = "PctEld",
              fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.greens"),
              fill.legend = tm_legend(title = "Persentase (%)"),
              col = "grey60", lwd = 0.2) +
  tm_title("Georgia: Persentase Penduduk Lansia (≥ 65 Tahun) (1990)") +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("left", "bottom")) +
  tm_layout(legend.outside = TRUE, frame = FALSE)
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Regresi Global (OLS) sebagai Titik Awal

Membangun Model OLS

# Model OLS (regresi global)
model_ols <- lm(PctBach ~ PctFB + PctPov + PctBlack + PctEld,
                data = georgia_sf)

# Ringkasan model
summary(model_ols)
## 
## Call:
## lm(formula = PctBach ~ PctFB + PctPov + PctBlack + PctEld, data = georgia_sf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1414  -2.0934  -0.2113   1.6709  19.9212 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.67113    1.63682   7.741 1.22e-12 ***
## PctFB        2.54521    0.29839   8.530 1.29e-14 ***
## PctPov      -0.28291    0.07810  -3.622 0.000396 ***
## PctBlack     0.07685    0.02803   2.742 0.006834 ** 
## PctEld      -0.10531    0.13784  -0.764 0.446047    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.013 on 154 degrees of freedom
## Multiple R-squared:  0.5163, Adjusted R-squared:  0.5037 
## F-statistic: 41.09 on 4 and 154 DF,  p-value: < 2.2e-16

Interpretasi output OLS:
Model OLS menghasilkan satu set koefisien yang berlaku untuk seluruh 159 kabupaten. Misalnya, koefisien PctPov menunjukkan rata-rata perubahan PctBach untuk setiap kenaikan 1% tingkat kemiskinan, diasumsikan sama di seluruh wilayah Georgia.

Diagnostik Regresi OLS

# Plot diagnostik residual
par(mfrow = c(2, 2))
plot(model_ols)

par(mfrow = c(1, 1))
# Uji normalitas residual (Shapiro-Wilk)
shapiro.test(residuals(model_ols))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_ols)
## W = 0.90536, p-value = 1.282e-08
# Uji heteroskedastisitas (Breusch-Pagan)
bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 40.41, df = 4, p-value = 3.56e-08
# Variance Inflation Factor (multikolinieritas)
vif(model_ols)
##    PctFB   PctPov PctBlack   PctEld 
## 1.334895 3.148163 2.328246 1.769013

Interpretasi VIF: - VIF < 5: Multikolinieritas dapat ditoleransi - VIF 5–10: Multikolinieritas moderat, perlu perhatian - VIF > 10: Multikolinieritas serius

Peta Residual OLS

# Simpan residual dan fitted value
georgia_map$resid_ols <- residuals(model_ols)
georgia_map$fitted_ols <- fitted(model_ols)

# Peta residual OLS
tm_shape(georgia_map) +
  tm_polygons(
    fill = "resid_ols",
    fill.scale = tm_scale_intervals(
      style = "jenks",
      n = 5,
      midpoint = 0,
      values = "brewer.rd_bu"
    ),
    fill.legend = tm_legend(title = "Residual OLS"),
    col = "grey60",
    col_alpha = 0.3,
    lwd = 0.2
  ) +
  tm_title("Peta Residual Model OLS") +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("left", "bottom")) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE
  )

Pertanyaan: Apakah residual tampak tersebar acak, atau ada pola spasial tertentu? Pola spasial dalam residual mengindikasikan bahwa OLS tidak sepenuhnya mampu menangkap efek spasial yang ada.

Uji Heterogenitas Spasial dan Autokorelasi Spasial

Uji Heterogenitas Spasial (Breusch-Pagan / Koenker-Bassett)

# Uji Breusch-Pagan (dari paket lmtest)
bp_test <- bptest(model_ols)
print(bp_test)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 40.41, df = 4, p-value = 3.56e-08
# Uji Koenker-Bassett yang lebih robust
# Menggunakan koordinat sebagai proksi spasial
georgia_coords <- st_coordinates(st_centroid(georgia_sf))
georgia_sf$x_coord <- georgia_coords[, 1]
georgia_sf$y_coord <- georgia_coords[, 2]

# Regresi residual kuadrat pada koordinat
resid_sq <- residuals(model_ols)^2
bp_spatial <- lm(resid_sq ~ x_coord + y_coord, data = georgia_sf)
summary(bp_spatial)
## 
## Call:
## lm(formula = resid_sq ~ x_coord + y_coord, data = georgia_sf)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -18.17 -14.70 -11.39  -4.09 379.59 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.387e+01  1.145e+02  -0.296    0.768
## x_coord      3.358e-06  3.652e-05   0.092    0.927
## y_coord      1.284e-05  2.776e-05   0.463    0.644
## 
## Residual standard error: 42.9 on 156 degrees of freedom
## Multiple R-squared:  0.001398,   Adjusted R-squared:  -0.0114 
## F-statistic: 0.1092 on 2 and 156 DF,  p-value: 0.8966

Mengapa OLS Tidak Cukup?

Bagian ini merangkum bukti empiris dari analisis di atas untuk memotivasi penggunaan GWR.

# Rangkuman hasil diagnostik OLS
cat("=== DIAGNOSTIK MODEL OLS ===\n\n")
## === DIAGNOSTIK MODEL OLS ===
cat("R-squared OLS:", round(summary(model_ols)$r.squared, 4), "\n")
## R-squared OLS: 0.5163
cat("   Adjusted R-squared:", round(summary(model_ols)$adj.r.squared, 4), "\n\n")
##    Adjusted R-squared: 0.5037
cat("Breusch-Pagan Test (Heteroskedastisitas):\n")
## Breusch-Pagan Test (Heteroskedastisitas):
cat("   p-value:", format.pval(bp_spatial$p.value, digits = 4), "\n\n")
##    p-value:

Pemilihan Fungsi Kernel

Pemilihan fungsi kernel mempengaruhi bagaimana bobot menurun seiring jarak. Tidak ada aturan baku — pemilihan bergantung pada konteks dan pengujian empiris.

Panduan Praktis Pemilihan Kernel

Kernel Karakteristik Kapan Digunakan
Bisquare Compact support, kurva lonceng terpotong Default yang direkomendasikan; efisien secara komputasi
Gaussian Semua pengamatan berkontribusi, penurunan halus Data sangat padat, tidak ada outlier spasial
Exponential Penurunan cepat dari pusat Efek lokal yang kuat, komponen jarak eksponensial
Tricube Compact support, lebih “rounded” dari bisquare Alternatif bisquare untuk permukaan yang lebih halus
Boxcar Semua pengamatan dalam radius = bobot sama Sangat jarang; untuk area yang homogen

Perbandingan Kernel Secara Empiris

# Perbandingan AICc untuk berbagai kernel (fixed bandwidth)
# Mengkonversi ke SpatialPolygonsDataFrame untuk GWmodel
georgia_sp <- as_Spatial(georgia_sf)

cat("Membandingkan AICc untuk berbagai fungsi kernel...\n")
## Membandingkan AICc untuk berbagai fungsi kernel...
cat("(Proses ini memerlukan beberapa menit)\n\n")
## (Proses ini memerlukan beberapa menit)
# Fungsi untuk menghitung GWR dengan kernel berbeda
compare_kernels <- function(kernel_type) {
  bw <- bw.gwr(PctBach ~ PctFB + PctPov + PctBlack + PctEld,
               data = georgia_sp,
               approach = "AICc",
               kernel = kernel_type,
               adaptive = TRUE,
               longlat = TRUE)
  
  gwr_fit <- gwr.basic(PctBach ~ PctFB + PctPov + PctBlack + PctEld,
                        data = georgia_sp,
                        bw = bw,
                        kernel = kernel_type,
                        adaptive = TRUE,
                        longlat = TRUE)
  
  return(list(kernel = kernel_type, bw = bw, 
              AICc = gwr_fit$GW.diagnostic$AICc,
              R2 = gwr_fit$GW.diagnostic$gw.R2))
}

# Jalankan perbandingan (bisa memakan waktu)
# kernels <- c("bisquare", "gaussian", "exponential", "tricube")
# hasil_kernel <- lapply(kernels, compare_kernels)
# 
# # Tampilkan hasil
# df_kernel <- data.frame(
#   Kernel = sapply(hasil_kernel, `[[`, "kernel"),
#   Bandwidth = sapply(hasil_kernel, `[[`, "bw"),
#   AICc = sapply(hasil_kernel, `[[`, "AICc"),
#   R2 = sapply(hasil_kernel, `[[`, "R2")
# )
# print(df_kernel[order(df_kernel$AICc), ])

# Untuk modul ini, kita gunakan bisquare (default yang direkomendasikan)
cat("Kernel yang dipilih: bisquare (compact support, direkomendasikan)\n")
## Kernel yang dipilih: bisquare (compact support, direkomendasikan)

Rekomendasi Praktis: Untuk sebagian besar kasus, gunakan kernel bisquare dengan adaptive bandwidth sebagai titik awal. Kernel ini efisien secara komputasi dan memberikan hasil yang stabil.

Pemilihan Bandwidth

Bandwidth adalah parameter krusial dalam GWR. Bandwidth yang terlalu kecil menghasilkan koefisien yang sangat fluktuatif (overfitting); bandwidth yang terlalu besar mendekati OLS global.

Metode Pemilihan Bandwidth

a. Cross-Validation (CV):
Meminimalkan prediksi leave-one-out cross-validation error: \[CV = \sum_{i=1}^{n} \left[y_i - \hat{y}_{\neq i}(h)\right]^2\]

b. AICc (direkomendasikan):
Meminimalkan Corrected Akaike Information Criterion yang menyeimbangkan kecocokan model dan kompleksitas.

Fixed vs. Adaptive Bandwidth

# ====================================================
# Pemilihan bandwidth: ADAPTIVE dengan AICc
# ====================================================
cat("Memilih bandwidth adaptif menggunakan AICc...\n")
## Memilih bandwidth adaptif menggunakan AICc...
bw_adaptive_AICc <- bw.gwr(
  PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data       = georgia_sp,
  approach   = "AICc",   # metode: "AICc" atau "CV"
  kernel     = "bisquare",
  adaptive   = TRUE,     # adaptive = jumlah tetangga
  longlat    = TRUE      # koordinat geografis (derajat)
)
## Adaptive bandwidth (number of nearest neighbours): 105 AICc value: 873.0522 
## Adaptive bandwidth (number of nearest neighbours): 73 AICc value: 872.1205 
## Adaptive bandwidth (number of nearest neighbours): 51 AICc value: 877.0335 
## Adaptive bandwidth (number of nearest neighbours): 84 AICc value: 871.3051 
## Adaptive bandwidth (number of nearest neighbours): 93 AICc value: 871.2177 
## Adaptive bandwidth (number of nearest neighbours): 97 AICc value: 871.4947 
## Adaptive bandwidth (number of nearest neighbours): 89 AICc value: 871.3171 
## Adaptive bandwidth (number of nearest neighbours): 94 AICc value: 871.4239 
## Adaptive bandwidth (number of nearest neighbours): 91 AICc value: 871.4916 
## Adaptive bandwidth (number of nearest neighbours): 93 AICc value: 871.2177
cat("Bandwidth adaptif optimal (AICc):", bw_adaptive_AICc, "tetangga\n\n")
## Bandwidth adaptif optimal (AICc): 93 tetangga
# ====================================================
# Pemilihan bandwidth: ADAPTIVE dengan CV
# ====================================================
cat("Memilih bandwidth adaptif menggunakan CV...\n")
## Memilih bandwidth adaptif menggunakan CV...
bw_adaptive_CV <- bw.gwr(
  PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data       = georgia_sp,
  approach   = "CV",
  kernel     = "bisquare",
  adaptive   = TRUE,
  longlat    = TRUE
)
## Adaptive bandwidth: 105 CV score: 2625.666 
## Adaptive bandwidth: 73 CV score: 2697.361 
## Adaptive bandwidth: 126 CV score: 2630.578 
## Adaptive bandwidth: 93 CV score: 2618.385 
## Adaptive bandwidth: 84 CV score: 2642.148 
## Adaptive bandwidth: 97 CV score: 2614.047 
## Adaptive bandwidth: 101 CV score: 2619.241 
## Adaptive bandwidth: 96 CV score: 2616.531 
## Adaptive bandwidth: 99 CV score: 2613.08 
## Adaptive bandwidth: 99 CV score: 2613.08
cat("Bandwidth adaptif optimal (CV):", bw_adaptive_CV, "tetangga\n\n")
## Bandwidth adaptif optimal (CV): 99 tetangga
# ====================================================
# Pemilihan bandwidth: FIXED dengan AICc
# ====================================================
cat("Memilih bandwidth fixed menggunakan AICc...\n")
## Memilih bandwidth fixed menggunakan AICc...
bw_fixed_AICc <- bw.gwr(
  PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data       = georgia_sp,
  approach   = "AICc",
  kernel     = "bisquare",
  adaptive   = FALSE,    # fixed = dalam satuan jarak
  longlat    = TRUE
)
## Fixed bandwidth: 12364.63 AICc value: 873.3223 
## Fixed bandwidth: 7643.289 AICc value: 878.4815 
## Fixed bandwidth: 15282.58 AICc value: 878.0422 
## Fixed bandwidth: 10561.24 AICc value: 871.6523 
## Fixed bandwidth: 9446.68 AICc value: 872.2683 
## Fixed bandwidth: 11250.07 AICc value: 871.8681 
## Fixed bandwidth: 10135.51 AICc value: 871.7524 
## Fixed bandwidth: 10824.35 AICc value: 871.6657 
## Fixed bandwidth: 10398.63 AICc value: 871.675 
## Fixed bandwidth: 10661.74 AICc value: 871.6497 
## Fixed bandwidth: 10723.85 AICc value: 871.6519 
## Fixed bandwidth: 10623.35 AICc value: 871.6497
cat("Bandwidth fixed optimal (AICc):", round(bw_fixed_AICc, 2), "km (approx)\n\n")
## Bandwidth fixed optimal (AICc): 10661.74 km (approx)
# Ringkasan perbandingan bandwidth
cat("=== RINGKASAN PEMILIHAN BANDWIDTH ===\n")
## === RINGKASAN PEMILIHAN BANDWIDTH ===
cat("Adaptive (AICc):", bw_adaptive_AICc, "tetangga\n")
## Adaptive (AICc): 93 tetangga
cat("Adaptive (CV)  :", bw_adaptive_CV, "tetangga\n")
## Adaptive (CV)  : 99 tetangga
cat("Fixed (AICc)   :", round(bw_fixed_AICc, 2), "km\n")
## Fixed (AICc)   : 10661.74 km
cat("\nPilihan untuk analisis: Adaptive bandwidth =", 
    bw_adaptive_AICc, "tetangga (AICc)\n")
## 
## Pilihan untuk analisis: Adaptive bandwidth = 93 tetangga (AICc)

Interpretasi Bandwidth Adaptif: Bandwidth = 50 tetangga berarti setiap estimasi lokal menggunakan 50 kabupaten terdekat. Dengan n = 159, ini menggunakan sekitar 31% data di setiap titik.

Visualisasi Profil Bandwidth

# Visualisasi perubahan AICc seiring bandwidth (jika didukung)
# Plot manual untuk memahami trade-off bandwidth
bw_range <- seq(30, 100, by = 5)

# Catatan: loop ini bisa memakan waktu; jalankan untuk n kecil
# aicc_vals <- sapply(bw_range, function(bw) {
#   gwr_fit <- gwr.basic(PctBach ~ PctFB + PctPov + PctBlack + PctEld,
#                         data = georgia_sp,
#                         bw = bw,
#                         kernel = "bisquare",
#                         adaptive = TRUE,
#                         longlat = TRUE)
#   return(gwr_fit$GW.diagnostic$AICc)
# })
# 
# plot(bw_range, aicc_vals, type = "b", pch = 19,
#      xlab = "Bandwidth (jumlah tetangga)",
#      ylab = "AICc",
#      main = "Profil AICc vs. Bandwidth")
# abline(v = bw_adaptive_AICc, col = "red", lty = 2)
# legend("topright", legend = "Bandwidth Optimal", 
#        col = "red", lty = 2)

Pemodelan GWR

Menjalankan Model GWR

# ====================================================
# MODEL GWR UTAMA
# ====================================================

gwr_model <- gwr.basic(
  formula  = PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data     = georgia_sp,
  bw       = bw_adaptive_AICc,
  kernel   = "bisquare",
  adaptive = TRUE,
  longlat  = TRUE,
  )

Melihat Output GWR

# Ringkasan komprehensif model GWR
print(gwr_model)
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2026-06-04 14:37:32.043969 
##    Call:
##    gwr.basic(formula = PctBach ~ PctFB + PctPov + PctBlack + PctEld, 
##     data = georgia_sp, bw = bw_adaptive_AICc, kernel = "bisquare", 
##     adaptive = TRUE, longlat = TRUE)
## 
##    Dependent (y) variable:  PctBach
##    Independent variables:  PctFB PctPov PctBlack PctEld
##    Number of data points: 159
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1414  -2.0934  -0.2113   1.6709  19.9212 
## 
##    Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
##    (Intercept) 12.67113    1.63682   7.741 1.22e-12 ***
##    PctFB        2.54521    0.29839   8.530 1.29e-14 ***
##    PctPov      -0.28291    0.07810  -3.622 0.000396 ***
##    PctBlack     0.07685    0.02803   2.742 0.006834 ** 
##    PctEld      -0.10531    0.13784  -0.764 0.446047    
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 4.013 on 154 degrees of freedom
##    Multiple R-squared: 0.5163
##    Adjusted R-squared: 0.5037 
##    F-statistic: 41.09 on 4 and 154 DF,  p-value: < 2.2e-16 
##    ***Extra Diagnostic information
##    Residual sum of squares: 2480.558
##    Sigma(hat): 3.974888
##    AIC:  900.0487
##    AICc:  900.6013
##    BIC:  789.8755
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 93 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Great Circle distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                   Min.   1st Qu.    Median   3rd Qu.    Max.
##    Intercept  5.615696  9.053463 11.135097 13.227625 17.6829
##    PctFB      1.359774  2.273029  2.684983  3.347688  4.7413
##    PctPov    -0.758644 -0.466808 -0.312448 -0.174473  0.2040
##    PctBlack  -0.042651  0.061736  0.084654  0.117802  0.1585
##    PctEld    -0.418820 -0.200138 -0.088642  0.298300  0.4124
##    ************************Diagnostic information*************************
##    Number of data points: 159 
##    Effective number of parameters (2trace(S) - trace(S'S)): 25.77023 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 133.2298 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 871.2177 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 842.7247 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 764.4913 
##    Residual sum of squares: 1646.393 
##    R-square value:  0.6789454 
##    Adjusted R-square value:  0.6163751 
## 
##    ***********************************************************************
##    Program stops at: 2026-06-04 14:37:32.049739

Output print(gwr_model) menampilkan beberapa bagian penting:

  1. GWR Formula: Formula yang digunakan
  2. Kernel dan Bandwidth: Jenis kernel dan nilai bandwidth
  3. Summary of GWR coefficient estimates: Statistik ringkasan (Min, Q1, Median, Q3, Max) untuk setiap koefisien lokal
  4. GW model diagnostic information: AICc, AIC, R² global, R² adj. global
  5. Results of model comparison: Perbandingan dengan OLS global
# Melihat komponen output
names(gwr_model)
## [1] "GW.arguments"  "GW.diagnostic" "lm"            "SDF"          
## [5] "timings"       "this.call"     "Ftests"
# Diagnostik model
gwr_model$GW.diagnostic
## $RSS.gw
## [1] 1646.393
## 
## $AIC
## [1] 842.7247
## 
## $AICc
## [1] 871.2177
## 
## $enp
## [1] 25.77023
## 
## $edf
## [1] 133.2298
## 
## $gw.R2
## [1] 0.6789454
## 
## $gwR2.adj
## [1] 0.6163751
## 
## $BIC
## [1] 764.4913

Ekstraksi Koefisien Lokal

# Ekstrak hasil GWR ke dataframe
gwr_results <- as.data.frame(gwr_model$SDF)
names(gwr_results)
##  [1] "Intercept"     "PctFB"         "PctPov"        "PctBlack"     
##  [5] "PctEld"        "y"             "yhat"          "residual"     
##  [9] "CV_Score"      "Stud_residual" "Intercept_SE"  "PctFB_SE"     
## [13] "PctPov_SE"     "PctBlack_SE"   "PctEld_SE"     "Intercept_TV" 
## [17] "PctFB_TV"      "PctPov_TV"     "PctBlack_TV"   "PctEld_TV"    
## [21] "Local_R2"      "coords.x1"     "coords.x2"
# Lihat beberapa baris pertama
head(gwr_results[, c("Intercept", "PctFB", "PctPov", 
                      "PctBlack", "PctEld", "Local_R2")], 10)
##    Intercept    PctFB     PctPov    PctBlack      PctEld  Local_R2
## 1  11.067358 2.142394 -0.6031003  0.14844580  0.37357210 0.7281565
## 2  12.952442 3.089092 -0.2657384  0.06784012 -0.15243942 0.5938488
## 3   8.790850 4.270273  0.1465487 -0.02533626 -0.34820654 0.7075223
## 4  14.947569 2.211147 -0.3827517  0.08099225 -0.10559219 0.5509759
## 5  15.999762 1.763246 -0.4417194  0.08850367 -0.09113261 0.5421244
## 6   9.668907 2.275880 -0.5545197  0.13756207  0.41241701 0.7137955
## 7  13.313892 2.905912 -0.2925062  0.07112972 -0.13168561 0.5776912
## 8  12.618980 3.247004 -0.2391053  0.06638177 -0.17938699 0.6069573
## 9  15.002758 1.694722 -0.5917670  0.13283858  0.07153385 0.6552784
## 10  9.370606 2.297304 -0.4993843  0.13615905  0.38510286 0.7235046
# Gabungkan hasil GWR dengan data spasial
georgia_gwr <- georgia_sf

# Koefisien lokal
georgia_gwr$coef_intercept <- gwr_results$Intercept
georgia_gwr$coef_PctFB     <- gwr_results$PctFB
georgia_gwr$coef_PctPov    <- gwr_results$PctPov
georgia_gwr$coef_PctBlack  <- gwr_results$PctBlack
georgia_gwr$coef_PctEld    <- gwr_results$PctEld

# R² lokal dan residual
georgia_gwr$local_R2       <- gwr_results$Local_R2
georgia_gwr$resid_gwr      <- gwr_results$residual
georgia_gwr$yhat_gwr       <- gwr_results$yhat

Statistik Ringkasan Koefisien Lokal

# Ringkasan distribusi koefisien lokal
cat("=== RINGKASAN KOEFISIEN LOKAL GWR ===\n\n")
## === RINGKASAN KOEFISIEN LOKAL GWR ===
coef_vars <- c("coef_intercept", "coef_PctFB", "coef_PctPov", 
               "coef_PctBlack", "coef_PctEld")

for (var in coef_vars) {
  vals <- st_drop_geometry(georgia_gwr)[[var]]
  cat(sprintf("%-20s | Min: %6.3f  Median: %6.3f  Max: %6.3f  SD: %6.3f\n",
              var, min(vals), median(vals), max(vals), sd(vals)))
}
## coef_intercept       | Min:  5.616  Median: 11.135  Max: 17.683  SD:  2.893
## coef_PctFB           | Min:  1.360  Median:  2.685  Max:  4.741  SD:  0.805
## coef_PctPov          | Min: -0.759  Median: -0.312  Max:  0.204  SD:  0.224
## coef_PctBlack        | Min: -0.043  Median:  0.085  Max:  0.158  SD:  0.048
## coef_PctEld          | Min: -0.419  Median: -0.089  Max:  0.412  SD:  0.259
cat("\n")
cat("Koefisien OLS global (pembanding):\n")
## Koefisien OLS global (pembanding):
print(round(coef(model_ols), 4))
## (Intercept)       PctFB      PctPov    PctBlack      PctEld 
##     12.6711      2.5452     -0.2829      0.0768     -0.1053

Evaluasi dan Perbandingan Model

Perbandingan Metrik Model

# Metrik OLS
r2_ols     <- summary(model_ols)$r.squared
r2adj_ols  <- summary(model_ols)$adj.r.squared
aic_ols    <- AIC(model_ols)
rss_ols    <- sum(residuals(model_ols)^2)

# Metrik GWR
r2_gwr     <- gwr_model$GW.diagnostic$gw.R2
r2adj_gwr  <- gwr_model$GW.diagnostic$gwR2.adj
aicc_gwr   <- gwr_model$GW.diagnostic$AICc
rss_gwr    <- sum(gwr_results$residual^2)

# Moran's I residual GWR
georgia_nb <- spdep::poly2nb(georgia_poly, queen = TRUE)
georgia_lw <- spdep::nb2listw(georgia_nb, style = "W", zero.policy = TRUE)
moran_gwr <- moran.test(georgia_gwr$resid_gwr, 
                         georgia_lw, 
                         zero.policy = TRUE)

#Moran's I residual OLS
moran_ols <- moran.test(residuals(model_ols), georgia_lw, zero.policy = TRUE)

cat("=== PERBANDINGAN OLS vs. GWR ===\n\n")
## === PERBANDINGAN OLS vs. GWR ===
cat(sprintf("%-30s %10s %10s\n", "Metrik", "OLS", "GWR"))
## Metrik                                OLS        GWR
cat(rep("-", 52), "\n", sep = "")
## ----------------------------------------------------
cat(sprintf("%-30s %10.4f %10.4f\n", "R-squared", r2_ols, r2_gwr))
## R-squared                          0.5163     0.6789
cat(sprintf("%-30s %10.4f %10.4f\n", "Adjusted R-squared", r2adj_ols, r2adj_gwr))
## Adjusted R-squared                 0.5037     0.6164
cat(sprintf("%-30s %10.2f %10.2f\n", "AIC / AICc", aic_ols, aicc_gwr))
## AIC / AICc                         900.05     871.22
cat(sprintf("%-30s %10.4f %10.4f\n", "RSS (Residual Sum of Sq)", rss_ols, rss_gwr))
## RSS (Residual Sum of Sq)        2480.5579  1646.3926
cat(sprintf("%-30s %10.4f %10.4f\n", "Moran's I Residual", 
            moran_ols$estimate[1], moran_gwr$estimate[1]))
## Moran's I Residual                 0.0378     0.0142
cat(sprintf("%-30s %10s %10s\n", "Moran's I p-value",
            ifelse(moran_ols$p.value < 0.001, "<0.001", round(moran_ols$p.value, 4)),
            ifelse(moran_gwr$p.value < 0.001, "<0.001", round(moran_gwr$p.value, 4))))
## Moran's I p-value                  0.1797     0.3356

Perbandingan Residual

# Visualisasi perbandingan residual OLS vs. GWR
# Visualisasi perbandingan residual OLS vs. GWR
par(mfrow = c(1, 2))

# PERBAIKAN: Mengubah georgia_sf menjadi georgia_map
hist(georgia_map$resid_ols, breaks = 20, 
     main = "Residual OLS", xlab = "Residual",
     col = "lightcoral", border = "white")
abline(v = 0, col = "red", lty = 2)

# Menggunakan objek gwr yang sudah sinkron
hist(georgia_gwr$resid_gwr, breaks = 20,
     main = "Residual GWR", xlab = "Residual",
     col = "lightblue", border = "white")
abline(v = 0, col = "blue", lty = 2)

par(mfrow = c(1, 1))
# Perbandingan fitted vs. actual
ggplot() +
  # OLS (PERBAIKAN: Menggunakan georgia_map)
  geom_point(data = st_drop_geometry(georgia_map),
             aes(x = PctBach, y = fitted_ols, color = "OLS"),
             alpha = 0.6, size = 2) +
  # GWR
  geom_point(data = st_drop_geometry(georgia_gwr),
             aes(x = PctBach, y = yhat_gwr, color = "GWR"),
             alpha = 0.6, size = 2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
  scale_color_manual(values = c("OLS" = "coral", "GWR" = "steelblue")) +
  labs(title = "Nilai Fitted vs. Aktual: OLS vs. GWR",
       x = "PctBach Aktual (%)",
       y = "PctBach Prediksi (%)",
       color = "Model") +
  theme_bw()

Interpretasi Koefisien Lokal

Prinsip Interpretasi

Koefisien lokal GWR \(\hat\beta_k(u_i, v_i)\) diinterpretasikan mirip dengan koefisien OLS, tetapi hanya berlaku untuk lokasi \((u_i, v_i)\):

“Di kabupaten X, setiap kenaikan 1% tingkat kemiskinan (PctPov) dikaitkan dengan perubahan sebesar \(\hat\beta_{PctPov}(u_X, v_X)\)% dalam persentase gelar sarjana (PctBach), dengan memperhitungkan variabel lain yang konstan.”

Tanda koefisien: - Positif (+): Kenaikan prediktor dikaitkan dengan kenaikan respons di lokasi tersebut - Negatif (−): Kenaikan prediktor dikaitkan dengan penurunan respons di lokasi tersebut - Mendekati nol: Variabel tidak berpengaruh signifikan secara lokal

Ringkasan Statistik Koefisien

# Tabel ringkasan koefisien lokal
coef_summary <- data.frame(
  Variabel = c("Intercept", "PctFB", "PctPov", "PctBlack", "PctEld"),
  OLS = round(coef(model_ols), 4),
  Min_GWR = round(apply(gwr_results[, c("Intercept","PctFB","PctPov",
                                          "PctBlack","PctEld")], 2, min), 4),
  Median_GWR = round(apply(gwr_results[, c("Intercept","PctFB","PctPov",
                                             "PctBlack","PctEld")], 2, median), 4),
  Max_GWR = round(apply(gwr_results[, c("Intercept","PctFB","PctPov",
                                          "PctBlack","PctEld")], 2, max), 4),
  SD_GWR = round(apply(gwr_results[, c("Intercept","PctFB","PctPov",
                                         "PctBlack","PctEld")], 2, sd), 4)
)

print(coef_summary)
##              Variabel     OLS Min_GWR Median_GWR Max_GWR SD_GWR
## (Intercept) Intercept 12.6711  5.6157    11.1351 17.6829 2.8930
## PctFB           PctFB  2.5452  1.3598     2.6850  4.7413 0.8054
## PctPov         PctPov -0.2829 -0.7586    -0.3124  0.2040 0.2242
## PctBlack     PctBlack  0.0768 -0.0427     0.0847  0.1585 0.0480
## PctEld         PctEld -0.1053 -0.4188    -0.0886  0.4124 0.2586
# Boxplot distribusi koefisien lokal
gwr_coef_long <- gwr_results %>%
  select(Intercept, PctFB, PctPov, PctBlack, PctEld) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Koefisien")

ggplot(gwr_coef_long, aes(x = Variabel, y = Koefisien, fill = Variabel)) +
  geom_boxplot(alpha = 0.7, outlier.size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  # Tambahkan garis untuk koefisien OLS
  geom_point(data = data.frame(
    Variabel = c("Intercept", "PctFB", "PctPov", "PctBlack", "PctEld"),
    OLS = coef(model_ols)
  ), aes(y = OLS), color = "black", shape = 18, size = 4) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Distribusi Koefisien Lokal GWR",
       subtitle = "Berlian hitam = koefisien OLS global",
       x = "Variabel Prediktor",
       y = "Nilai Koefisien Lokal") +
  theme_bw() +
  theme(legend.position = "none")

Pemetaan Koefisien Lokal

Peta Koefisien Lokal untuk Setiap Prediktor

# 1. Aktifkan mode plot statis
tmap_mode("plot")

# ====================================================
# SINKRONISASI DATA
# Memasukkan kolom koefisien dari georgia_gwr ke georgia_map
# ====================================================
georgia_map$coef_PctFB    <- georgia_gwr$coef_PctFB
georgia_map$coef_PctPov   <- georgia_gwr$coef_PctPov
georgia_map$coef_PctBlack <- georgia_gwr$coef_PctBlack
georgia_map$coef_PctEld   <- georgia_gwr$coef_PctEld

# ====================================================
# PETA CHOROPLETH 1: % Foreign Born
# ====================================================
map_coef_PctFB <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "coef_PctFB",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(
      title = "β(PctFB)", 
      position = tm_pos_in("left", "bottom"), 
      frame = TRUE
    )
  ) +
  tm_title("Koefisien GWR: % Foreign Born", size = 0.8)

# ====================================================
# PETA CHOROPLETH 2: % Kemiskinan
# ====================================================
map_coef_PctPov <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "coef_PctPov",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(
      title = "β(PctPov)", 
      position = tm_pos_in("left", "bottom"), 
      frame = TRUE
    )
  ) +
  tm_title("Koefisien GWR: % Kemiskinan", size = 0.8)

# ====================================================
# PETA CHOROPLETH 3: % Black
# ====================================================
map_coef_PctBlack <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "coef_PctBlack",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(
      title = "β(PctBlack)", 
      position = tm_pos_in("left", "bottom"), 
      frame = TRUE
    )
  ) +
  tm_title("Koefisien GWR: % Black", size = 0.8)

# ====================================================
# PETA CHOROPLETH 4: % Lansia
# ====================================================
map_coef_PctEld <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "coef_PctEld",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(
      title = "β(PctEld)", 
      position = tm_pos_in("left", "bottom"), 
      frame = TRUE
    )
  ) +
  tm_title("Koefisien GWR: % Lansia", size = 0.8)

# ====================================================
# TAMPILKAN 4 PETA CHOROPLETH BERDAMPINGAN
# ====================================================
tmap_arrange(map_coef_PctFB, map_coef_PctPov, 
             map_coef_PctBlack, map_coef_PctEld,
             ncol = 2)

Berdasarkan analisis, kita dapat memberikan interpretasi seperti berikut:

Contoh interpretasi (bergantung pada hasil aktual):

  • PctPov (% Kemiskinan): Secara global (OLS), setiap kenaikan 1% kemiskinan dikaitkan dengan penurunan ~X% gelar sarjana. Namun, GWR menunjukkan bahwa di kabupaten-kabupaten metropolitan (Atlanta dan sekitarnya), pengaruh kemiskinan lebih kuat negatif dibanding di wilayah perdesaan.

  • PctFB (% Foreign Born): Koefisien lokal bervariasi dari negatif di bagian tertentu hingga positif di bagian lain, menunjukkan bahwa pengaruh imigran terhadap pendidikan tidak seragam.

  • PctBlack (% Black): Pengaruh komposisi etnis bervariasi, mencerminkan heterogenitas historis dalam akses pendidikan antar-wilayah.

Peta R² Lokal

georgia_map$local_R2 <- georgia_gwr$local_R2

# ====================================================
# Peta R² Lokal (Choropleth / Poligon Murni tmap v4)
# ====================================================
map_local_R2 <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "local_R2",                      # Menggunakan fill untuk poligon
    fill.scale = tm_scale_intervals(        # Mengelompokkan parameter skala
      style = "jenks",
      n = 5,
      values = "brewer.yl_gn_bu"            # Menggunakan standar nama palet v4
    ),
    fill.legend = tm_legend(
      title = "R² Lokal",
      position = tm_pos_in("left", "bottom"), # Legenda di dalam agar aman dari bug layout
      frame = TRUE
    )
  ) +
  tm_title("GWR: R² Lokal (Kecocokan Model per Kabupaten)", size = 0.8)

# Tampilkan peta
map_local_R2

Interpretasi R² Lokal: Area dengan R² lokal tinggi menunjukkan bahwa variabel prediktor yang dipilih mampu menjelaskan variasi PctBach dengan baik di area tersebut. R² lokal rendah mengindikasikan bahwa ada variabel penting yang tidak termasuk dalam model, atau hubungannya sangat tidak linear.

Peta Residual GWR

georgia_map$resid_gwr <- georgia_gwr$resid_gwr

# ====================================================
# Peta Residual GWR (Choropleth / Poligon Murni tmap v4)
# ====================================================
map_resid_gwr <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "resid_gwr",                     # Mewarnai area poligon berdasarkan nilai residual
    fill.scale = tm_scale_intervals(        # Manajemen skala interval di tmap v4
      style = "jenks",
      n = 5,
      values = "brewer.rd_bu",              # Menggunakan standar nama palet v4
      midpoint = 0                          # Pusat warna (putih/netral) dikunci di angka 0
    ),
    fill.legend = tm_legend(
      title = "Residual GWR",
      position = tm_pos_in("left", "bottom"), # Legenda di dalam area peta agar aman dari bug layout
      frame = TRUE
    )
  ) +
  tm_title("Peta Residual Model GWR", size = 0.8)

# Tampilkan peta
map_resid_gwr

13.4 Peta Perbandingan: Sebelum dan Sesudah GWR

georgia_map$resid_ols <- as.numeric(residuals(model_ols))
georgia_map$resid_gwr <- as.numeric(georgia_gwr$resid_gwr)

semua_residual <- c(georgia_map$resid_ols, georgia_map$resid_gwr)
breaks_bersama <- pretty(semua_residual, n = 5)

# Pembuatan peta perbandingan multivariate facet (tmap v4)
peta_perbandingan <- tm_shape(georgia_map) +
  tm_polygons(
    fill = c("resid_ols", "resid_gwr"),
    fill.scale = tm_scale_intervals(
      breaks = breaks_bersama, 
      values = "brewer.rd_bu", 
      midpoint = 0
    ),
    fill.legend = tm_legend(
      title = "Nilai Residual",
      position = tm_pos_out("right", "center")
    )
  ) +
  tm_facets(ncol = 2)

peta_perbandingan
## [plot mode] legend/component: Some components or legends are too "high" and are
## therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Uji Signifikansi Koefisien Lokal

Pengantar Uji Signifikansi Lokal

Koefisien lokal GWR disertai dengan standard error dan t-value, tetapi interpretasi langsung bermasalah karena: 1. Multiple testing problem: Kita melakukan uji hipotesis di 159 lokasi 2. Dependensi spasial: Pengujian di lokasi yang berdekatan tidak independen

Paket GWmodel menyediakan fungsi gwr.t.adjust() yang mengkoreksi masalah uji berganda menggunakan metode Byrne et al. (2016).

Ekstrak Standard Error dan T-values Lokal

# Ekstrak standard error dan t-values dari model GWR
# (diperlukan se.fit = TRUE saat menjalankan gwr.basic)

# Nama kolom standard error
se_cols <- paste0(c("Intercept", "PctFB", "PctPov", "PctBlack", "PctEld"), "_SE")
tv_cols <- paste0("TV_", c("Intercept", "PctFB", "PctPov", "PctBlack", "PctEld"))

# Cek nama kolom yang tersedia
names(gwr_results)
##  [1] "Intercept"     "PctFB"         "PctPov"        "PctBlack"     
##  [5] "PctEld"        "y"             "yhat"          "residual"     
##  [9] "CV_Score"      "Stud_residual" "Intercept_SE"  "PctFB_SE"     
## [13] "PctPov_SE"     "PctBlack_SE"   "PctEld_SE"     "Intercept_TV" 
## [17] "PctFB_TV"      "PctPov_TV"     "PctBlack_TV"   "PctEld_TV"    
## [21] "Local_R2"      "coords.x1"     "coords.x2"
# T-values lokal
if ("PctPov_TV" %in% names(gwr_results)) {
  georgia_gwr$tv_PctFB    <- gwr_results$PctFB_TV
  georgia_gwr$tv_PctPov   <- gwr_results$PctPov_TV
  georgia_gwr$tv_PctBlack <- gwr_results$PctBlack_TV
  georgia_gwr$tv_PctEld   <- gwr_results$PctEld_TV
  
  cat("T-values lokal berhasil diekstrak.\n")
}
## T-values lokal berhasil diekstrak.

Koreksi Uji Berganda: Metode gwr.t.adjust

# Koreksi p-value dengan multiple testing adjustment
gwr_sig <- gwr.t.adjust(gwr_model)

# Lihat struktur output
summary(gwr_sig)
##         Length Class                  Mode
## results   6    -none-                 list
## SDF     159    SpatialPointsDataFrame S4
# Ekstrak p-values terkoreksi ke data spasial
# Metode Bonferroni atau FDR (tergantung versi GWmodel)

# Menambahkan hasil signifikansi ke dataframe
sig_results <- as.data.frame(gwr_sig$SDF)

if ("PctPov_p" %in% names(sig_results) || 
    "PctPov_pval" %in% names(sig_results)) {
  
  # Cek nama kolom p-value
  p_col_check <- names(sig_results)[grepl("PctPov", names(sig_results))]
  cat("Kolom p-value tersedia:", p_col_check, "\n")
}
## Kolom p-value tersedia: PctPov_t PctPov_p PctPov_p_by PctPov_p_fb PctPov_p_bo PctPov_p_bh

Klasifikasi Signifikansi

# Kategori signifikansi berdasarkan t-value
# Aturan umum: |t| > 1.96 → p < 0.05
# Untuk multiple testing: gunakan threshold lebih ketat

categorize_sig <- function(t_val) {
  case_when(
    abs(t_val) > 2.576 ~ "Sangat signifikan (p<0.01)",
    abs(t_val) > 1.960 ~ "Signifikan (p<0.05)",
    abs(t_val) > 1.645 ~ "Marjinal (p<0.10)",
    TRUE               ~ "Tidak signifikan"
  )
}

if ("PctPov_TV" %in% names(gwr_results)) {
  georgia_gwr$sig_PctPov   <- categorize_sig(gwr_results$PctPov_TV)
  georgia_gwr$sig_PctFB    <- categorize_sig(gwr_results$PctFB_TV)
  georgia_gwr$sig_PctBlack <- categorize_sig(gwr_results$PctBlack_TV)
  georgia_gwr$sig_PctEld   <- categorize_sig(gwr_results$PctEld_TV)
  
  # Ringkasan signifikansi
  cat("\n=== RINGKASAN SIGNIFIKANSI KOEFISIEN LOKAL ===\n\n")
  for (var in c("sig_PctFB", "sig_PctPov", "sig_PctBlack", "sig_PctEld")) {
    cat(var, ":\n")
    print(table(st_drop_geometry(georgia_gwr)[[var]]))
    cat("\n")
  }
}
## 
## === RINGKASAN SIGNIFIKANSI KOEFISIEN LOKAL ===
## 
## sig_PctFB :
## 
## Sangat signifikan (p<0.01) 
##                        159 
## 
## sig_PctPov :
## 
##          Marjinal (p<0.10) Sangat signifikan (p<0.01) 
##                          6                         90 
##        Signifikan (p<0.05)           Tidak signifikan 
##                         21                         42 
## 
## sig_PctBlack :
## 
##          Marjinal (p<0.10) Sangat signifikan (p<0.01) 
##                         24                         54 
##        Signifikan (p<0.05)           Tidak signifikan 
##                         29                         52 
## 
## sig_PctEld :
## 
##   Marjinal (p<0.10) Signifikan (p<0.05)    Tidak signifikan 
##                  33                  18                 108

Pemetaan Signifikansi

Peta Koefisien dengan Masking Signifikansi

Pendekatan umum adalah menampilkan koefisien lokal dengan stippling (titik-titik) atau masking untuk menunjukkan mana yang signifikan dan mana yang tidak.

# Sinkronisasi data ke poligon & filter signifikansi (|t| > 1.96)
georgia_map$tv_PctPov <- gwr_results$PctPov_TV
georgia_map$significant <- abs(georgia_map$tv_PctPov) > 1.96

# Pembuatan peta t-value choropleth (tmap v4)
map_tvalue <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "tv_PctPov",
    fill.scale = tm_scale_intervals(style = "jenks", n = 7, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(title = "t-value (PctPov)", position = tm_pos_out("right"))
  ) +
  tm_shape(georgia_map[georgia_map$significant, ]) +
  tm_borders(col = "black", lwd = 1.8) +
  tm_title("T-values Koefisien Lokal: % Kemiskinan\n(Batas tebal = p < 0.05)", size = 0.85)

map_tvalue

# Membuat kategori gabungan arah + signifikansi (menggunakan PctPov_TV dari hasil sebelumnya)
# Membuat 3 kategori signifikansi berdasarkan t-value
georgia_map$PctPov_sig_cat3 <- case_when(
  gwr_results$PctPov_TV < -1.96 ~ "Negatif Signifikan",
  gwr_results$PctPov_TV >  1.96 ~ "Positif Signifikan",
  TRUE                          ~ "Tidak Signifikan"
)

# Pembuatan peta choropleth 3 kategori (tmap v4)
map_sig_cat3 <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "PctPov_sig_cat3",
    fill.scale = tm_scale_categorical(values = c(
      "Negatif Signifikan" = "#4575b4",  # Biru
      "Positif Signifikan" = "#d73027",  # Merah
      "Tidak Signifikan"   = "#e0e0e0"   # Abu-abu netral
    )),
    fill.legend = tm_legend(title = "Signifikansi PctPov", position = tm_pos_out("right"))
  ) +
  tm_title("Peta Koefisien × Signifikansi: % Kemiskinan", size = 0.9)

map_sig_cat3

Peta R² Lokal dengan Koefisien

# Sinkronisasi data ke objek poligon georgia_map
georgia_map$local_R2 <- georgia_gwr$local_R2
georgia_map$coef_PctPov <- georgia_gwr$coef_PctPov

# Peta 1: R² Lokal
map1 <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "local_R2",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.yl_gn_bu"),
    fill.legend = tm_legend(title = "R² Lokal", position = tm_pos_out("right"))
  ) +
  tm_title("R² Lokal GWR", size = 0.85)

# Peta 2: Koefisien PctPov
map2 <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "coef_PctPov",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.rd_bu", midpoint = 0),
    fill.legend = tm_legend(title = "β(PctPov)", position = tm_pos_out("right"))
  ) +
  tm_title("Koefisien Lokal: % Kemiskinan", size = 0.85)

# Tampilkan kedua peta berdampingan dengan ukuran sama besar
tmap_arrange(map1, map2, ncol = 2)

Keterbatasan GWR

Meskipun GWR adalah alat yang sangat berguna, penting untuk memahami keterbatasannya:

Multikolinieritas Lokal

Kolinieritas merupakan masalah yang dapat terjadi pada model regresi mana pun, namun urgensinya menjadi jauh lebih tinggi dalam konteks GWR. Hal ini disebabkan karena model GWR mengestimasi ulang parameter pada subset geografis lokal. Pada skala lokal tersebut, peubah-peubah penjelas berpotensi besar memiliki korelasi yang jauh lebih kuat dibandingkan pada skala wilayah studi secara keseluruhan (mengikuti prinsip spasial bahwa hal-hal yang mirip cenderung saling berdekatan). Multikolinieritas lokal ini dapat memicu inflasi varians (variance inflation), yang berdampak pada ketidakstabilan model dan menghasilkan estimasi parameter yang tidak andal atau tidak presisi.

Sebagai langkah awal, kita melakukan diagnostik kolinieritas pada hasil model. Jika pengujian menunjukkan adanya bukti multikolinieritas yang mengkhawatirkan, kita memiliki dua pilihan solusi: mengubah spesifikasi model dengan menghapus peubah penjelas penentu (culprit covariates), atau beralih ke bentuk model yang kebal (robust) terhadap kolinieritas, seperti local ridge regression yang tersedia melalui fungsi gwr.lcr(). Untuk saat ini, kita akan fokus pada tahap diagnostik terlebih dahulu menggunakan fungsi gwr.collin.diagno() yang menghasilkan Condition Number (CN) lokal serta nilai Variance Inflation Factor (VIF) lokal dari model yang telah fit.

# 1. Jalankan fungsi diagnostik kolinieritas lokal GWR
local_cn <- gwr.collin.diagno(
  formula = PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data = georgia_sp,
  bw = 93,
  kernel = "bisquare",
  adaptive = TRUE,
  longlat = TRUE
)

# 2. Tampilkan ringkasan statistik Condition Number lokal di Console
summary(local_cn$SDF$local_CN)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.40   13.93   15.26   15.68   16.87   23.51
# 3. Sinkronisasi nilai CN ke data poligon asli sebagai vektor numerik
georgia_map$cond_num <- as.numeric(local_cn$SDF$local_CN)

# 4. Visualisasi peta choropleth Condition Number menggunakan tmap v4
map_cn <- tm_shape(georgia_map) +
  tm_polygons(
    fill = "cond_num",
    fill.scale = tm_scale_intervals(style = "jenks", n = 5, values = "brewer.or_rd"),
    fill.legend = tm_legend(title = "Condition Number", position = tm_pos_out("right"))
  ) +
  tm_title("Multikolinieritas Lokal: Condition Number", size = 0.85)

# Tampilkan peta
map_cn

Interpretasi OutputEvaluasi Nilai Ambang Kritis (Threshold): + Berdasarkan nilai \(\eta_i\), jika \(\eta_i > 30\), hal tersebut mengindikasikan adanya multikolinieritas lokal tingkat moderat. Jika \(\eta_i > 100\), matriks lokal mendekati singular (terjadi multikolinieritas lokal yang sangat serius).

  • Analisis Nilai Maksimum (Summary): Perhatikan nilai Max pada luaran summary(local_cn\(SDF\)local_CN). Jika nilai maksimum tersebut berada di bawah angka 30, maka pembagi \(\lambda_{\min}\) pada seluruh lokasi tidak ada yang mendekati nol, menandakan struktur data lokal di seluruh kabupaten Georgia aman dari gangguan kolinieritas.

  • Pola Distribusi Spasial (Peta): Melalui peta choropleth, wilayah dengan warna jingga tua hingga merah menunjukkan area di mana peubah penjelas saling bertumpang tindih secara spasial. Jika wilayah dengan nilai CN > 30 muncul, estimasi koefisien lokal \(\beta(u_i, v_i)\) di area tersebut menjadi tidak stabil, sehingga disarankan untuk menggunakan pendekatan local ridge regression (gwr.lcr()) khusus pada kluster wilayah tersebut.

# ====================================================
# 1. CETAK VIF GLOBAL (OLS)
# ====================================================
cat("=== DIAGNOSTIK MULTIKOLINIERITAS GLOBAL (OLS) ===\n")
## === DIAGNOSTIK MULTIKOLINIERITAS GLOBAL (OLS) ===
vif_global <- car::vif(model_ols)
print(round(vif_global, 3))
##    PctFB   PctPov PctBlack   PctEld 
##    1.335    3.148    2.328    1.769
cat("\n")
# ====================================================
# 2. CETAK CONDITION NUMBER LOKAL (GWR)
# ====================================================
cat("=== DIAGNOSTIK MULTIKOLINIERITAS LOKAL (GWR) ===\n")
## === DIAGNOSTIK MULTIKOLINIERITAS LOKAL (GWR) ===
summary(local_cn$SDF$local_CN)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.40   13.93   15.26   15.68   16.87   23.51

VIF Global (OLS): Mengukur tingkat keterikatan linear antara satu peubah penjelas tertentu (\(X_k\)) dengan peubah-peubah penjelas lainnya di dalam model. Hasilnya adalah satu nilai untuk setiap peubah \(X\), yang berlaku sama untuk seluruh wilayah studi.

Condition Number Lokal (GWR): Mengukur tingkat ketidakstabilan atau kepelikan numerik dari seluruh sistem matriks koordinat di lokasi spesifik \(i\). Hasilnya adalah satu nilai untuk setiap titik lokasi/kabupaten (\(i\)), yang merangkum kondisi seluruh peubah \(X\) setelah dikenakan pembobotan geografis.

Ringkasan Keterbatasan

Keterbatasan Penjelasan Solusi/Mitigasi
Multikolinieritas lokal Variabel lebih berkorelasi di sub-wilayah Periksa condition number lokal; hapus variabel kolinear
Multiple testing Uji hipotesis di banyak lokasi sekaligus Gunakan gwr.t.adjust() untuk koreksi
Overfitting Bandwidth kecil → model terlalu menyesuaikan noise Pilih bandwidth via AICc; validasi silang
Outlier sensitif Pengamatan ekstrem sangat memengaruhi estimasi lokal Diagnostik leverage; robust GWR
Interpretasi Koefisien bervariasi ≠ hubungan kausal berbeda GWR bersifat eksploratif, bukan kausal
Ekstrapolasi Koefisien di tepi wilayah kurang dapat diandalkan Gunakan lebih banyak data di tepi; catat ketidakpastian
Bandwidth stationarity Satu bandwidth untuk semua variabel Pertimbangkan Mixed GWR atau MGWR

Mixed GWR (MGWR): Pengembangan GWR yang memungkinkan bandwidth berbeda untuk setiap variabel prediktor. Jika beberapa variabel bersifat global dan lainnya lokal, Mixed GWR lebih tepat. Tersedia melalui gwr.mixed() di GWmodel.

Latihan Mandiri

Latihan 1: Variasi Kernel

Jalankan model GWR dengan tiga jenis kernel berbeda (bisquare, gaussian, tricube) menggunakan bandwidth yang sama. Bandingkan AICc dan pola spasial koefisien yang dihasilkan. Kernel mana yang menghasilkan AICc terkecil?

Latihan 2: Fixed vs. Adaptive Bandwidth

Jalankan GWR dengan bandwidth fixed (hasil dari bw_fixed_AICc) dan bandingkan hasilnya dengan bandwidth adaptive. Jelaskan perbedaan dalam pola koefisien lokal yang dihasilkan.

Latihan 3: Variabel Respons Berbeda

Ulangi seluruh analisis menggunakan MedInc (median pendapatan) sebagai variabel respons dan PctBach, PctFB, PctPov, PctBlack sebagai prediktor. Apakah pola non-stasioneritas sama?

Latihan 4: Mixed GWR

Coba jalankan Mixed GWR menggunakan gwr.mixed() dengan asumsi bahwa PctEld adalah variabel global dan variabel lainnya lokal:

# Mixed GWR (beberapa variabel global, lainnya lokal)
gwr_mixed <- gwr.mixed(
  PctBach ~ PctFB + PctPov + PctBlack + PctEld,
  data          = georgia_sp,
  fixed.vars    = "PctEld",  # variabel global
  intercept.fixed = FALSE,
  bw            = bw_adaptive_AICc,
  kernel        = "bisquare",
  adaptive      = TRUE,
  longlat       = TRUE
)
print(gwr_mixed)

Data Anda Sendiri

Cari dataset spasial dari wilayah Indonesia (misalnya: data Badan Pusat Statistik tingkat kecamatan atau kabupaten) dan terapkan seluruh alur kerja GWR. Beberapa sumber data: - BPS: https://www.bps.go.id/ - Satu Data Indonesia: https://data.go.id/ - PODES (Potensi Desa): data sosial-ekonomi tingkat desa

Referensi

Comber, L. (2019). INIA International Workshop on Spatial Analysis in R — Session 7: Regression and Spatial Heterogeneity Effects with Geographically Weighted Regression (GWR). RPubs. https://rpubs.com/lexcomber/INIA_Session7

GWmodel Team (2016). Geographically Weighted Regression Practical. RPubs. https://rpubs.com/gwmodel/176883

Kram, M. (2023). EPI 563: Spatial Epidemiology — Week 12: Spatial Regression III: Geographically Weighted Regression. https://mkram01.github.io/EPI563-SpatialEPI/spatial-regression-iii-geographically-weighted-regression.html

Potier, C. (2018). Lab 3: Geographically Weighted Regression — A Case Study of Children’s Social Skill Scores in Vancouver, Canada. GEOB479, University of British Columbia.

Bivand, R. (2022). Geographically Weighted Regression (Package Vignette: spgwr). CRAN. https://cran.r-project.org/web/packages/spgwr/vignettes/GWR.pdf