Model GWR

1. Pengenalan Geographically Weighted Regression (GWR)

Geographically Weighted Regression (GWR) adalah metode regresi yang memperhitungkan variasi spasial dalam data. GWR memungkinkan parameter model regresi berubah di setiap lokasi, memberikan fleksibilitas untuk menangkap pola lokal daripada satu model global yang seragam.

1.1 Kapan Menggunakan GWR?

GWR cocok digunakan ketika terdapat indikasi bahwa hubungan antara variabel dependen dan independen berubah-ubah berdasarkan lokasi geografis. GWR sering digunakan dalam penelitian yang melibatkan data spasial seperti ekologi, geografi, epidemiologi, dan ekonomi regional.

2. Model GWR

Model GWR dituliskan sebagai:

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

Di mana: - \(y_i\) adalah nilai variabel dependen di lokasi \(i\). - \(\beta_k(u_i, v_i)\) adalah parameter regresi untuk variabel \(k\) pada lokasi geografis \((u_i, v_i)\). - \(x_{ik}\) adalah variabel independen ke-\(k\) pada lokasi \(i\). - \(\varepsilon_i\) adalah error residual pada lokasi \(i\).

Setiap lokasi \((u_i, v_i)\) memiliki nilai \(\beta_k\) yang berbeda, membuat GWR model lokal di setiap titik.

3. Pemilihan Bandwidth dalam GWR

Bandwidth dalam GWR menentukan seberapa banyak pengaruh lokasi tetangga di sekitar titik yang diamati. Pilihan bandwidth sangat mempengaruhi hasil estimasi parameter. Terdapat dua jenis bandwidth yang umum digunakan:

  1. Fixed Bandwidth: Semua titik observasi menggunakan ukuran bandwidth yang sama, baik untuk lokasi yang padat maupun jarang.

  2. Adaptive Bandwidth: Bandwidth disesuaikan dengan kepadatan titik observasi, sehingga di wilayah dengan observasi lebih sedikit, bandwidth diperluas untuk memasukkan lebih banyak tetangga.

3.1 Metode Pemilihan Bandwidth

Pemilihan bandwidth yang optimal dilakukan dengan metode berikut:

  • Cross-Validation (CV): Bandwidth dipilih yang meminimalkan nilai fungsi loss pada data yang tidak dilatih.

  • Akaike Information Criterion (AICc): Bandwidth dipilih berdasarkan nilai AIC yang disesuaikan, dengan memilih bandwidth yang meminimalkan AIC.

4. Metode Estimasi dalam GWR

Untuk mengestimasi parameter dalam model GWR, digunakan metode weighted least squares (WLS). Setiap titik observasi memiliki bobot yang berbeda berdasarkan jaraknya dari titik target. Secara matematis, estimasi \(\beta_k(u_i, v_i)\) diestimasi dengan:

\[ \hat{\beta}(u_i, v_i) = (X^T W_i X)^{-1} X^T W_i y \]

Di mana: - \(W_i\) adalah matriks bobot untuk lokasi ke-\(i\), dan elemen dari matriks ini ditentukan oleh fungsi kernel yang bergantung pada jarak antara titik \(i\) dan titik lainnya. - \(X\) adalah matriks desain (covariates), dan \(y\) adalah vektor variabel dependen.

4.1 Fungsi Kernel

Bobot dihitung menggunakan fungsi kernel yang mendekati nol ketika jarak semakin besar. Beberapa jenis kernel yang umum digunakan: - Gaussian Kernel: \[ w_{ij} = \exp \left( - \frac{d_{ij}^2}{2h^2} \right) \] - Bisquare Kernel: \[ w_{ij} = \left( 1 - \left( \frac{d_{ij}}{h} \right)^2 \right)^2 \text{ jika } d_{ij} < h, \text{ sebaliknya } w_{ij} = 0 \]

5. Pengujian dalam GWR

5.1 Pengujian Signifikansi Global

Untuk menguji apakah GWR lebih baik daripada model regresi global, dapat digunakan test statistik F. Pengujian ini membandingkan varians dari model global dengan varians yang dihasilkan oleh GWR.

5.2 Pengujian Signifikansi Lokal

Di setiap lokasi, koefisien \(\beta_k(u_i, v_i)\) dapat diuji secara lokal. t-statistik dapat dihitung untuk setiap koefisien pada setiap lokasi, guna melihat apakah parameter tersebut signifikan secara statistik pada level tertentu.

6. Teori Uji Asumsi dalam GWR

Sama seperti regresi linear biasa, GWR memiliki asumsi-asumsi yang harus diuji. Namun, beberapa asumsi ini memiliki nuansa berbeda karena adanya faktor spasial.

6.1 Uji Normalitas Residual

Normalitas residual dapat diuji dengan metode Shapiro-Wilk Test. Hipotesis nol dari uji ini adalah bahwa residual berdistribusi normal.

Rumus uji Shapiro-Wilk:

\[ W = \frac{\left( \sum_{i=1}^{n} a_i r_{(i)} \right)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]

Di mana: - \(r_{(i)}\) adalah residual yang diurutkan. - \(a_i\) adalah koefisien yang diperoleh dari matriks kovarian teoretis distribusi normal. - \(\bar{y}\) adalah rata-rata residual.

Jika nilai \(p\) kecil (misalnya, kurang dari 0.05), maka dapat disimpulkan bahwa residual tidak berdistribusi normal.

6.2 Uji Homoskedastisitas

Untuk memeriksa homoskedastisitas (atau varians yang tetap), Breusch-Pagan Test dapat digunakan. Uji ini memeriksa apakah error residual bervariasi secara sistematis terhadap variabel prediktor.

Rumus uji Breusch-Pagan:

\[ BP = \frac{n}{2} R^2 \]

Di mana: - \(R^2\) adalah koefisien determinasi dari model regresi tambahan yang mengkuadratkan residual sebagai variabel dependen. - \(n\) adalah jumlah observasi.

Hipotesis nol menyatakan bahwa varians residual konstan (homoskedastisitas). Jika nilai \(p\) kecil, ada indikasi heteroskedastisitas.

6.3 Uji Autokorelasi Spasial

Uji autokorelasi spasial, seperti Moran’s I, digunakan untuk memeriksa apakah residual memiliki pola spasial yang signifikan.

Rumus Moran’s I:

\[ I = \frac{n}{W} \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij} (y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]

Di mana: - \(n\) adalah jumlah observasi. - \(w_{ij}\) adalah elemen dari matriks bobot spasial. - \(\bar{y}\) adalah rata-rata nilai residual.

Nilai \(I\) yang mendekati +1 menunjukkan autokorelasi positif, sementara nilai mendekati -1 menunjukkan autokorelasi negatif. Jika \(I\) signifikan, berarti ada pola spasial dalam residual.

6.4 Uji Multikolinearitas

Untuk mendeteksi multikolinearitas antara variabel independen, dapat digunakan Variance Inflation Factor (VIF).

Rumus VIF:

\[ VIF_k = \frac{1}{1 - R_k^2} \]

Di mana: - \(R_k^2\) adalah koefisien determinasi dari regresi variabel \(X_k\) terhadap variabel independen lainnya.

Jika VIF lebih dari 10, hal ini menunjukkan adanya multikolinearitas yang serius, sehingga model perlu dikoreksi.

7. Kesimpulan

GWR adalah alat yang kuat untuk analisis data spasial karena mampu menangkap variasi lokal dalam hubungan antara variabel. Dengan pemilihan bandwidth yang tepat, estimasi yang akurat, dan pengujian asumsi yang lengkap, GWR dapat memberikan wawasan yang lebih dalam tentang bagaimana hubungan antara variabel berubah di ruang geografis.

APLIKASI - R

Data Generating

Derajat (Latitude longitude)

data_spasial <- data.frame(
  y = runif(100),
  x1 = runif(100),
  x2 = runif(100),
  longitude = runif(100, min = -100, max = -90),
  latitude = runif(100, min = 30, max = 40)
)

coordinates(data_spasial) <- ~longitude+latitude

Optimal Bandwidth

bandwidth_optimal <- gwr.sel(
  y ~ x1 + x2,
  data = data_spasial,
  coords = cbind(data_spasial$longitude, data_spasial$latitude),
  adapt = TRUE
)
## Warning in gwr.sel(y ~ x1 + x2, data = data_spasial, coords =
## cbind(data_spasial$longitude, : data is Spatial* object, ignoring coords
## argument
## Adaptive q: 0.381966 CV score: 9.77187 
## Adaptive q: 0.618034 CV score: 9.660247 
## Adaptive q: 0.763932 CV score: 9.614546 
## Adaptive q: 0.854102 CV score: 9.596684 
## Adaptive q: 0.9098301 CV score: 9.580846 
## Adaptive q: 0.9442719 CV score: 9.575263 
## Adaptive q: 0.9869169 CV score: 9.566139 
## Adaptive q: 0.970628 CV score: 9.569923 
## Adaptive q: 0.9806951 CV score: 9.567951 
## Adaptive q: 0.9919142 CV score: 9.564653 
## Adaptive q: 0.9950027 CV score: 9.563681 
## Adaptive q: 0.9969115 CV score: 9.563096 
## Adaptive q: 0.9980912 CV score: 9.562741 
## Adaptive q: 0.9988203 CV score: 9.562524 
## Adaptive q: 0.9992709 CV score: 9.562391 
## Adaptive q: 0.9995494 CV score: 9.562309 
## Adaptive q: 0.9997215 CV score: 9.562258 
## Adaptive q: 0.9998279 CV score: 9.562227 
## Adaptive q: 0.9998936 CV score: 9.562207 
## Adaptive q: 0.9999343 CV score: 9.562195 
## Adaptive q: 0.9999343 CV score: 9.562195
## Warning in gwr.sel(y ~ x1 + x2, data = data_spasial, coords =
## cbind(data_spasial$longitude, : Bandwidth converged to upper bound:1

Running GWR

model_gwr <- gwr(
  y ~ x1 + x2,
  data = data_spasial,
  coords = cbind(data_spasial$longitude, data_spasial$latitude),
  adapt = bandwidth_optimal,
  hatmatrix = TRUE,
  se.fit = TRUE
)
## Warning in gwr(y ~ x1 + x2, data = data_spasial, coords =
## cbind(data_spasial$longitude, : data is Spatial* object, ignoring coords
## argument

Koefisien pada setiap lokasi

coef_gwr <- model_gwr$SDF@data[, c("x1", "x2")]

summary

# Menampilkan ringkasan model
summary(model_gwr)
##           Length Class                  Mode     
## SDF         100  SpatialPointsDataFrame S4       
## lhat      10000  -none-                 numeric  
## lm           11  -none-                 list     
## results      14  -none-                 list     
## bandwidth   100  -none-                 numeric  
## adapt         1  -none-                 numeric  
## hatmatrix     1  -none-                 logical  
## gweight       1  -none-                 character
## gTSS          1  -none-                 numeric  
## this.call     7  -none-                 call     
## fp.given      1  -none-                 logical  
## timings      12  -none-                 numeric

UTM (Universal Transverse Mercator)

Mendefinisikan CRS latitude/longitude (WGS84)

proj_longlat <- CRS("+proj=longlat +datum=WGS84")

Mendefinisikan CRS UTM, misalnya UTM zona 33N

proj_utm <- CRS("+proj=utm +zone=33 +datum=WGS84")

Definisikan CRS untuk data longitude/latitude (WGS84)

proj4string(data_spasial) <- CRS("+proj=longlat +datum=WGS84")

Konversi data spatial ke UTM

data_spasial_utm <- spTransform(data_spasial, proj_utm)

Menjalankan GWR dengan data yang sudah terproyeksi ke UTM

bandwidth_optimal <- gwr.sel(y ~ x1 + x2, data = data_spasial_utm, adapt = TRUE)
## Adaptive q: 0.381966 CV score: 9.795876 
## Adaptive q: 0.618034 CV score: 9.666732 
## Adaptive q: 0.763932 CV score: 9.622607 
## Adaptive q: 0.9271036 CV score: 9.584867 
## Adaptive q: 0.8647776 CV score: 9.599351 
## Adaptive q: 0.9549476 CV score: 9.576879 
## Adaptive q: 0.9721561 CV score: 9.573872 
## Adaptive q: 0.9827915 CV score: 9.56956 
## Adaptive q: 0.9893646 CV score: 9.568099 
## Adaptive q: 0.9965049 CV score: 9.565563 
## Adaptive q: 0.9937775 CV score: 9.566544 
## Adaptive q: 0.9978399 CV score: 9.565095 
## Adaptive q: 0.998665 CV score: 9.564809 
## Adaptive q: 0.9991749 CV score: 9.564634 
## Adaptive q: 0.9994901 CV score: 9.564526 
## Adaptive q: 0.9996848 CV score: 9.56446 
## Adaptive q: 0.9998052 CV score: 9.564419 
## Adaptive q: 0.9998796 CV score: 9.564394 
## Adaptive q: 0.9999256 CV score: 9.564378 
## Adaptive q: 0.9999256 CV score: 9.564378
## Warning in gwr.sel(y ~ x1 + x2, data = data_spasial_utm, adapt = TRUE):
## Bandwidth converged to upper bound:1
model_gwr <- gwr(y ~ x1 + x2, data = data_spasial_utm, adapt = bandwidth_optimal)

UJI ASUMSI

Heterogeneity

Model OLS untuk perbandingan

model_ols <- lm(y ~ x1 + x2, data = data_spasial)
summary(model_ols)
## 
## Call:
## lm(formula = y ~ x1 + x2, data = data_spasial)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49645 -0.27224  0.00619  0.24107  0.58561 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.41301    0.08632   4.785 6.12e-06 ***
## x1           0.14568    0.09585   1.520    0.132    
## x2          -0.01803    0.11607  -0.155    0.877    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.304 on 97 degrees of freedom
## Multiple R-squared:  0.02476,    Adjusted R-squared:  0.004653 
## F-statistic: 1.231 on 2 and 97 DF,  p-value: 0.2964
bptest(model_ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 0.44035, df = 2, p-value = 0.8024

Lakukan perbandingan AIC (Semakin rendah AIC, semakin baik model)

Jika AIC model GWR lebih rendah dari model OLS, maka ada heterogenitas spasial.

AIC(model_ols)
## [1] 50.56841
model_gwr$results$AICc
## NULL

Multikolinearity

Gunakan VIF untuk mendeteksi multikolinearitas

vif(model_ols)
##       x1       x2 
## 1.023671 1.023671

Uji Spatial Autokorelasi

Ekstraksi residual GWR

residual_gwr <- model_gwr$SDF$gwr.e

Membuat matriks bobot spasial berdasarkan tetangga terdekat

coords <- cbind(data_spasial$longitude, data_spasial$latitude)
nb <- knn2nb(knearneigh(coords, k=4))
listw <- nb2listw(nb, style="W")

Uji Moran’s I untuk residual spasial

moran.test(residual_gwr, listw)
## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw    
## 
## Moran I statistic standard deviate = -0.30648, p-value = 0.6204
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.030064622      -0.010101010       0.004242966

Uji Normality

shapiro.test(residual_gwr)
## 
##  Shapiro-Wilk normality test
## 
## data:  residual_gwr
## W = 0.95571, p-value = 0.002012

Membuat plot -spplot

Ekstrak koefisien lokal untuk variabel x1 dan x2

coef_x1 <- model_gwr$SDF$x1
coef_x2 <- model_gwr$SDF$x2

Plot koefisien regresi lokal untuk variabel x1

spplot(model_gwr$SDF, "x1", main = "Koefisien GWR untuk x1", col.regions = terrain.colors(100))

Plot koefisien regresi lokal untuk variabel x2

spplot(model_gwr$SDF, "x2", main = "Koefisien GWR untuk x2", col.regions = terrain.colors(100))

Membuat plot -gglot2

Konversi SpatialPointsDataFrame menjadi data frame biasa

gwr_df <- as.data.frame(model_gwr$SDF)

# Plot koefisien lokal untuk x1
ggplot(gwr_df, aes(x = longitude, y = latitude)) +
  geom_point(aes(color = x1), size = 2) +
  scale_color_viridis_c() +
  labs(title = "Koefisien GWR untuk x1", color = "Koefisien") +
  theme_minimal()

# Plot koefisien lokal untuk x2
ggplot(gwr_df, aes(x = longitude, y = latitude)) +
  geom_point(aes(color = x2), size = 2) +
  scale_color_viridis_c() +
  labs(title = "Koefisien GWR untuk x2", color = "Koefisien") +
  theme_minimal()

Membuat Surface Plot

IDW

Assuming data_spasial is in WGS84

If data_spasial is in UTM, you may want to transform it to WGS84

proj4string(data_spasial)
## [1] "+proj=longlat +datum=WGS84 +no_defs"
data_spasial_wgs84 <- spTransform(data_spasial, CRS("+proj=longlat +datum=WGS84 +no_defs"))

Create the grid with the same CRS

Tentukan batas (extent) dari peta

x.range <- range(data_spasial$longitude)
y.range <- range(data_spasial$latitude)

Buat grid (titik-titik prediksi di mana kita ingin mengestimasi nilai koefisien)

grd <- expand.grid(
  longitude = seq(from = x.range[1], to = x.range[2], length.out = 100),
  latitude = seq(from = y.range[1], to = y.range[2], length.out = 100)
)

coordinates(grd) <- ~longitude + latitude
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

Gabungkan koefisien dengan data spasial

data_spasial$coef_x1 <- coef_x1

Lakukan IDW untuk variabel koefisien x1

idw_output <- idw(coef_x1 ~ 1, data_spasial, newdata = grd, idp = 2.0)  # idp = power parameter
## [inverse distance weighted interpolation]

Konversi output IDW menjadi data frame untuk ggplot

idw_df <- as.data.frame(idw_output)
ggplot() +
  geom_tile(data = idw_df, aes(x = longitude, y = latitude, fill = var1.pred)) +
  scale_fill_viridis_c() +
  labs(title = "Surface Plot of GWR Coefficients for x1 (IDW)",
       fill = "Koefisien x1") +
  theme_minimal()