UTS Spasial 2024
Silvi
2024-10-17
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.
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.
Load the necessary libraries
library(spdep) ## Warning: package 'spdep' was built under R version 4.2.3## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Warning: package 'sf' was built under R version 4.2.3## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(rgdal) ## Warning: package 'rgdal' was built under R version 4.2.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.2.3## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/LENOVO/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/LENOVO/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.1-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(raster) ## Warning: package 'raster' was built under R version 4.2.3library(corrplot) ## Warning: package 'corrplot' was built under R version 4.2.3## corrplot 0.92 loadedlibrary(DescTools) ## Warning: package 'DescTools' was built under R version 4.2.3library(nortest)
library(car) ## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:DescTools':
##
## Recodelibrary(spatialreg)## Warning: package 'spatialreg' was built under R version 4.2.3## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.2.3##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(sf)
library(lmtest)## Warning: package 'lmtest' was built under R version 4.2.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'stringr' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.0
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ dplyr::select() masks raster::select()
## ✖ purrr::some() masks car::some()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(spgwr)## Warning: package 'spgwr' was built under R version 4.2.3## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)library(GWmodel)## Warning: package 'GWmodel' was built under R version 4.2.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.2.3## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 4.2.3## Welcome to GWmodel version 2.3-2.library(dplyr)DESCRIPTIVE ANALYSIS
1.1 Set Working Directory
setwd("D:/kuliah/semester 5/Spatial")1.2 Load Data
data <- read.csv("Data Spatial UTS New.csv", sep=";", dec=",")
jabar <- read_sf("D:/kuliah/semester 5/Spatial/jabar_modified.shp")
# Convert shapefile to spatial object
CoordK <- as_Spatial(jabar)
# Check summary of the data
summary(data)## Kabupaten.Kota y x1 x2
## Length:27 Min. : 924 Min. : 382.0 Min. : 33661
## Class :character 1st Qu.: 3804 1st Qu.: 824.5 1st Qu.: 92713
## Mode :character Median : 5384 Median : 1449.0 Median :216758
## Mean : 7527 Mean : 3873.3 Mean :238810
## 3rd Qu.:10414 3rd Qu.: 5749.0 3rd Qu.:324124
## Max. :26912 Max. :15047.0 Max. :701147
## x3 longitude latitude
## Min. : 2647 Min. :-7701397 Min. : 1081971
## 1st Qu.: 14858 1st Qu.:-6994542 1st Qu.:106878230
## Median : 22198 Median :-6836531 Median :107512302
## Mean : 29535 Mean :-6605407 Mean : 92565328
## 3rd Qu.: 41564 3rd Qu.:-6515690 3rd Qu.:108016295
## Max. :110195 Max. : -687121 Max. :108557818#####MAPPING####
jabar$id<-c(1:27)
data$id<-c(1:27)
jabar_sf<-st_as_sf(jabar)
jabar_merged <- jabar_sf %>%
left_join(data, by = "id")# Menghitung centroid dari setiap kabupaten
jabar_centroid <- st_centroid(jabar_merged)## Warning: st_centroid assumes attributes are constant over geometriesggplot() +
geom_sf(data = jabar_sf, aes(fill = y), color = "white") +
geom_sf(data = jabar_centroid, color = "red", size = 2) + # Menambahkan titik centroid sebagai warna merah
theme_bw() +
scale_fill_gradient(low = "pink", high = "brown") +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right",
axis.text.x = element_blank(),
axis.text.y = element_blank()
)
1.3 Explolatory Data
1.3.1 Boxplot variabel y (Penyakit TBC)
boxplot(data$y, main="Sebaran Peubah Y")
1.3.2 Boxplot variabel X1 (Kepadatan Penduduk)
boxplot(data$x1, main="Sebaran Peubah X1")
1.3.3 Boxplot variabel X2 (RT ber PHBS)
boxplot(data$x2, main="Sebaran Peubah X2")
1.3.4 Boxplot variabel X3 (Balita Imunisasi BCG)
boxplot(data$x3, main="Sebaran Peubah X3")
SPATIAL ANALYSIS
2.1 Spatial distribution plot for each variable
2.1.1 Spatial distribution plot for variable Y
k <- 27
colfunc <- colorRampPalette(c("green","yellow","red"))
color <- colfunc(k)
# Spatial distribution of Y
CoordK$y <- data$y
spplot(CoordK, "y", col.regions=color, main="Sebaran Spasial Peubah Y")
# Categorize data$y based on breaks and labels, with intervals including the upper bound
jabar_merged$y <- cut(data$y,
breaks = c(-Inf, 3804, 5384, 10414, Inf),
labels = c("Very Low", "Low", "High", "Very High"),
right = TRUE)
# Plot with green-yellow-red color scheme
ggplot() +
geom_sf(data = jabar_merged, aes(fill = y), color = NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "#00FF00", # Green
"Low" = "#FFFF00", # Yellow
"High" = "#FFA500", # Orange
"Very High" = "#FF0000" # Red
)) +
labs(title = "Tuberkulosis", fill = "y") +
theme(legend.position = "right",
axis.text.x = element_blank(),
axis.text.y = element_blank())
2.1.2 Spatial distribution plot for variabel X1
# Spatial distribution of X1
CoordK$x1 <- data$x1
spplot(CoordK, "x1", col.regions=color, main="Sebaran Spasial Peubah X1")
2.1.3 Spatial distribution plot for variabel X2
# Spatial distribution of X2
CoordK$x2 <- data$x2
spplot(CoordK, "x2", col.regions=color, main="Sebaran Spasial Peubah X2")
2.1.4 Spatial distribution plot for variabel X3
# Spatial distribution of X3
CoordK$x3 <- data$x3
spplot(CoordK, "x3", col.regions=color, main="Sebaran Spasial Peubah X3")
Linear Regression (OLS)
reg.klasik <- lm(y ~ x1 + x2 + x3, data=data)
summary(reg.klasik)##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4156.3 -1351.9 -159.7 1067.7 5080.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.106e+03 1.074e+03 1.029 0.3141
## x1 2.761e-01 1.205e-01 2.292 0.0314 *
## x2 -8.052e-03 4.949e-03 -1.627 0.1173
## x3 2.463e-01 3.330e-02 7.397 1.6e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2457 on 23 degrees of freedom
## Multiple R-squared: 0.8338, Adjusted R-squared: 0.8122
## F-statistic: 38.47 on 3 and 23 DF, p-value: 3.949e-093.1 Assumption tests for OLS Multicollinearity check using VIF
car::vif(reg.klasik)## x1 x2 x3
## 1.337374 3.304786 2.8550283.2 Normality test for residuals
shapiro.test(reg.klasik$residuals)##
## Shapiro-Wilk normality test
##
## data: reg.klasik$residuals
## W = 0.97304, p-value = 0.68343.3 Heteroskedasticity test (Breusch-Pagan)
bptest(reg.klasik)##
## studentized Breusch-Pagan test
##
## data: reg.klasik
## BP = 12.395, df = 3, p-value = 0.0061453.4 Asumsi AUtokorelasi
# Install dan load paket yang diperlukan
install.packages("spdep")## Warning: package 'spdep' is in use and will not be installedlibrary(spdep)
coords <- coordinates(CoordK) # koord untuk titik atau poligon
nb <- knn2nb(knearneigh(coords, k = 4)) # contoh dengan k = 4 tetangga terdekat
listw <- nb2listw(nb, style = "W") # Membuat matriks bobot spasial
# Menghitung Moran's I
# Misalkan variabel spasial yang dianalisis adalah 'y' dalam objek CoordK
moran_test <- moran.test(CoordK$y, listw)
# Menampilkan hasil uji Moran's I
print(moran_test)##
## Moran I test under randomisation
##
## data: CoordK$y
## weights: listw
##
## Moran I statistic standard deviate = 2.4656, p-value = 0.00684
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.23568399 -0.03846154 0.01236306Spatial Weigths Matrices
4.1 Calculate spatial weights matrices (W)
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_k}(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.2 Spatial Weight Matrices K-NN with k=5
longlat <- coordinates(CoordK)
W.knn <- knn2nb(knearneigh(longlat, k=5, longlat=TRUE))
W.knn1 <- nb2listw(W.knn, style="W")
# Plot KNN with k=5
plot(CoordK, col="pink", border="red", main="KNN dengan k=5")
plot(W.knn1, coords=longlat, col="blue", lwd=2, add=TRUE)
4.3 Spatial Weight Matrices Queen
W.queen <- nb2listw(poly2nb(CoordK, queen=TRUE), style="W", zero.policy=TRUE)
MI.queen <- moran.test(data$y, W.queen, randomisation=TRUE, zero.policy=TRUE)4.4 Test spatial dependence using Moran’s I for different spatial weight matrices
# Display Moran's I
MI.queen$estimate## Moran I statistic Expectation Variance
## 0.37947921 -0.03846154 0.01659271longlat <- coordinates(CoordK)
W.queen <- nb2listw(poly2nb(CoordK, queen=TRUE), style="W", zero.policy=TRUE)
W.knn <- knn2nb(knearneigh(longlat, k=5, longlat=TRUE))
W.knn1 <- nb2listw(W.knn, style="W")MORAN’S I
5.1 Calculate Global Moran’s I for both weights
MI.global.queen <- moran.test(data$y, W.queen, randomisation=TRUE, zero.policy=TRUE)
MI.global.knn <- moran.test(data$y, W.knn1, randomisation=TRUE, zero.policy=TRUE)5.2 Global Moran’s I Queen
MI.global.queen$estimate## Moran I statistic Expectation Variance
## 0.37947921 -0.03846154 0.016592715.3 Global Moran’s I K-NN
MI.global.knn$estimate## Moran I statistic Expectation Variance
## 0.268590785 -0.038461538 0.0095014345.4 Compare and select the best
if (MI.global.queen$estimate[1] > MI.global.knn$estimate[1]) {
best_weight <- "Queen"
best_moran <- MI.global.queen$estimate[1]
} else {
best_weight <- "KNN"
best_moran <- MI.global.knn$estimate[1]
}
cat("Matriks bobot terbaik adalah:", best_weight, "dengan nilai Moran's I global:", best_moran, "\n")## Matriks bobot terbaik adalah: Queen dengan nilai Moran's I global: 0.3794792LOCAL MORAN’S I
6.1 Calculate Local Moran’s I using the best weight
if (best_weight == "Queen") {
local_moran <- localmoran(data$y, W.queen, zero.policy=TRUE)
} else {
local_moran <- localmoran(data$y, W.knn1, zero.policy=TRUE)
}6.2 Create a data frame for local Moran’s I results
local_moran_df <- data.frame(
id = 1:nrow(data),
Local_Moran_I = local_moran[, 1],
p_value = local_moran[, 5]
)# Filter significant results
significant_local_moran <- local_moran_df %>% filter(p_value < 0.05)
# Display significant results
print(significant_local_moran)## id Local_Moran_I p_value
## 1 1 1.541908696 0.0097884291
## 7 7 0.673806488 0.0210302073
## 16 16 1.216855420 0.0010089682
## 19 19 1.609735205 0.0005572489
## 23 23 1.502639125 0.0096735085
## 24 24 -0.003240921 0.0012100779FUNGSI KERNEL
Bobot dihitung menggunakan fungsi kernel yang mendekati nol ketika jarak semakin besar. Beberapa jenis kernel yang umum digunakan:
7.1 Geographically Weighted Regression (GWR) Bi-square kernel
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 \]
b.gwr <- gwr.sel(y ~ x1 + x2 + x3, data=data, coords=longlat, gweight=gwr.bisquare)## Bandwidth: 0.894596 CV score: 340328471
## Bandwidth: 1.446042 CV score: 241023300
## Bandwidth: 1.786854 CV score: 233820597
## Bandwidth: 1.675766 CV score: 235188159
## Bandwidth: 1.891569 CV score: 233211953
## Bandwidth: 1.935724 CV score: 233257389
## Bandwidth: 1.902451 CV score: 233208900
## Bandwidth: 1.900574 CV score: 233208696
## Bandwidth: 1.900256 CV score: 233208691
## Bandwidth: 1.900215 CV score: 233208691
## Bandwidth: 1.900175 CV score: 233208691
## Bandwidth: 1.900215 CV score: 233208691gwr.model <- gwr(y ~ x1 + x2 + x3, data=data, coords=longlat, gweight=gwr.bisquare, bandwidth=b.gwr, hatmatrix=TRUE)
gwr.model## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data, coords = longlat,
## bandwidth = b.gwr, gweight = gwr.bisquare, hatmatrix = TRUE)
## Kernel function: gwr.bisquare
## Fixed bandwidth: 1.900215
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -3.2438e+02 2.7146e+01 9.7214e+02 2.4062e+03 3.6795e+03
## x1 5.7055e-02 1.7137e-01 3.1296e-01 4.6238e-01 5.1469e-01
## x2 -1.5095e-02 -1.2030e-02 -7.5229e-03 -1.3997e-03 7.8870e-04
## x3 1.8974e-01 2.0926e-01 2.4433e-01 2.6065e-01 2.6883e-01
## Global
## X.Intercept. 1105.5696
## x1 0.2761
## x2 -0.0081
## x3 0.2463
## Number of data points: 27
## Effective number of parameters (residual: 2traceS - traceS'S): 7.601862
## Effective degrees of freedom (residual: 2traceS - traceS'S): 19.39814
## Sigma (residual: 2traceS - traceS'S): 2311.78
## Effective number of parameters (model: traceS): 6.29636
## Effective degrees of freedom (model: traceS): 20.70364
## Sigma (model: traceS): 2237.706
## Sigma (ML): 1959.497
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 507.0322
## AIC (GWR p. 96, eq. 4.22): 492.263
## Residual sum of squares: 103669964
## Quasi-global R2: 0.8758941gwr.model$SDF@data## sum.w (Intercept) x1 x2 x3 (Intercept)_se
## 1 12.97326 3588.98803 0.06603811 -0.0146046204 0.2676420 1516.333
## 2 11.78972 3043.35203 0.11547311 -0.0124897015 0.2628005 1548.746
## 3 16.18408 1570.11770 0.24710507 -0.0092223544 0.2543056 1111.639
## 4 18.96005 809.60880 0.32917551 -0.0067978699 0.2426597 1000.871
## 5 17.02599 275.72741 0.39373242 -0.0042995028 0.2316739 1031.994
## 6 14.53501 -220.08696 0.47381138 -0.0009571322 0.2094606 1151.954
## 7 14.24373 -244.79898 0.49041569 -0.0001171626 0.1987284 1198.324
## 8 13.90377 -159.53460 0.48314951 -0.0003770142 0.1972992 1187.966
## 9 13.60218 27.16640 0.46233481 -0.0014037536 0.2004614 1139.482
## 10 16.85266 306.70087 0.41267767 -0.0038540352 0.2186882 1044.096
## 11 18.89810 619.12774 0.36248291 -0.0058840391 0.2329801 1004.694
## 12 15.55205 735.76008 0.36208385 -0.0061944880 0.2281014 1036.410
## 13 18.04816 1291.62353 0.28397313 -0.0089568345 0.2466971 1042.095
## 14 18.39747 1677.88219 0.23965895 -0.0101637202 0.2539881 1084.700
## 15 15.23304 2315.96867 0.17980595 -0.0123415079 0.2601415 1217.288
## 16 13.68218 3000.42925 0.11643539 -0.0140393084 0.2658436 1352.496
## 17 18.94065 1349.48794 0.27077465 -0.0088645440 0.2508612 1039.299
## 18 10.97654 -323.36454 0.51469496 0.0007886955 0.1897404 1319.423
## 19 13.37694 3431.15129 0.08009455 -0.0142458063 0.2668638 1484.003
## 20 14.67815 2496.35650 0.16293013 -0.0117180331 0.2611603 1311.974
## 21 19.62044 972.14484 0.31296334 -0.0075228897 0.2443257 1002.561
## 22 13.56491 27.12531 0.46243515 -0.0013956542 0.2003639 1140.113
## 23 13.27174 3337.83605 0.08650969 -0.0146409409 0.2677177 1422.827
## 24 12.73099 3679.49910 0.05705484 -0.0150949683 0.2688296 1508.258
## 25 19.46636 1147.12898 0.29305610 -0.0081952361 0.2475082 1015.149
## 26 15.36016 -153.83561 0.46806103 -0.0012239015 0.2090549 1134.272
## 27 12.55918 -324.37774 0.50569928 0.0006479104 0.1925244 1278.454
## x1_se x2_se x3_se gwr.e pred pred.se localR2
## 1 0.1354692 0.005281296 0.03349361 1518.59613 25393.4039 1831.0055 0.8445412
## 2 0.1364069 0.005430082 0.03368367 2132.16780 8600.8322 737.5168 0.8479551
## 3 0.1169601 0.004814673 0.03182909 2746.51623 6233.4838 1233.3016 0.8534678
## 4 0.1138947 0.004638513 0.03125113 -1478.09041 13773.0904 774.3149 0.8568215
## 5 0.1191180 0.004966019 0.03384499 -343.11788 8744.1179 690.2534 0.8579169
## 6 0.1322235 0.005967040 0.04489775 632.84967 2842.1503 780.4108 0.8543923
## 7 0.1386655 0.006615450 0.05167225 -401.33606 3332.3361 770.9978 0.8511516
## 8 0.1406492 0.006839704 0.05342069 328.45104 3321.5490 633.7242 0.8450568
## 9 0.1383037 0.006654160 0.05189902 -1357.66369 9285.6637 1115.7282 0.8366805
## 10 0.1241737 0.005383316 0.03936283 498.56765 3687.4323 593.6130 0.8459357
## 11 0.1163499 0.004784140 0.03311391 -258.02397 3195.0240 639.4692 0.8517679
## 12 0.1199720 0.005122589 0.03609793 -1477.74129 6511.7413 1069.6630 0.8402510
## 13 0.1135895 0.004607238 0.03093501 858.06558 4691.9344 587.5131 0.8476930
## 14 0.1145245 0.004615423 0.03084749 1023.36732 4360.6327 646.7027 0.8492262
## 15 0.1194676 0.004808662 0.03152309 3656.96941 8712.0306 598.3186 0.8412114
## 16 0.1262908 0.005030094 0.03253237 -2848.54396 14107.5440 1401.8186 0.8396613
## 17 0.1135652 0.004600823 0.03088515 -2213.93206 6172.9321 526.5483 0.8537553
## 18 0.1472024 0.007485853 0.05619418 -43.85476 967.8548 917.6991 0.8730387
## 19 0.1335973 0.005232378 0.03331385 3405.93272 6690.0673 822.9215 0.8446406
## 20 0.1249991 0.005040610 0.03261824 -1267.82751 4672.8275 842.4086 0.8482880
## 21 0.1131163 0.004578557 0.03092304 4669.71261 12827.2874 1192.7975 0.8549703
## 22 0.1384131 0.006665189 0.05199194 -745.65569 4846.6557 996.0625 0.8365475
## 23 0.1303233 0.005143008 0.03303674 -1832.95641 15394.9564 1305.1208 0.8407857
## 24 0.1352668 0.005274653 0.03353089 -2793.73862 10312.7386 1142.4786 0.8427819
## 25 0.1129842 0.004570217 0.03078139 -2549.79530 6977.7953 1148.1565 0.8543590
## 26 0.1315132 0.005923354 0.04468881 730.88554 3840.1145 662.3424 0.8525802
## 27 0.1448800 0.007259116 0.05520811 123.90228 1016.0977 984.2165 0.8620642
## (Intercept)_se_EDF x1_se_EDF x2_se_EDF x3_se_EDF pred.se
## 1 1566.528 0.1399535 0.005456119 0.03460233 1891.6161
## 2 1600.013 0.1409223 0.005609831 0.03479868 761.9303
## 3 1148.437 0.1208317 0.004974050 0.03288271 1274.1268
## 4 1034.002 0.1176649 0.004792059 0.03228562 799.9465
## 5 1066.155 0.1230611 0.005130406 0.03496534 713.1023
## 6 1190.086 0.1366004 0.006164563 0.04638396 806.2442
## 7 1237.991 0.1432556 0.006834437 0.05338273 796.5196
## 8 1227.290 0.1453050 0.007066114 0.05518903 654.7020
## 9 1177.202 0.1428819 0.006874428 0.05361700 1152.6614
## 10 1078.658 0.1282842 0.005561517 0.04066583 613.2630
## 11 1037.952 0.1202013 0.004942506 0.03421005 660.6371
## 12 1070.718 0.1239433 0.005292159 0.03729285 1105.0713
## 13 1076.591 0.1173496 0.004759748 0.03195903 606.9611
## 14 1120.606 0.1183155 0.004768204 0.03186862 668.1100
## 15 1257.583 0.1234222 0.004967840 0.03256658 618.1243
## 16 1397.266 0.1304713 0.005196602 0.03360926 1448.2221
## 17 1073.702 0.1173245 0.004753121 0.03190752 543.9782
## 18 1363.098 0.1520752 0.007733652 0.05805434 948.0771
## 19 1533.127 0.1380197 0.005405582 0.03441661 850.1621
## 20 1355.403 0.1291369 0.005207466 0.03369798 870.2943
## 21 1035.748 0.1168607 0.004730118 0.03194666 1232.2819
## 22 1177.853 0.1429949 0.006885823 0.05371299 1029.0344
## 23 1469.926 0.1346373 0.005313253 0.03413034 1348.3233
## 24 1558.185 0.1397445 0.005449256 0.03464084 1180.2973
## 25 1048.753 0.1167243 0.004721501 0.03180032 1186.1631
## 26 1171.819 0.1358666 0.006119430 0.04616811 684.2674
## 27 1320.774 0.1496759 0.007499409 0.05703563 1016.7964Plot beta coefficients for GWR model
CoordK@data$b1 <- gwr.model$SDF$x1
spplot(CoordK, zcol="b1", main="Peta Sebaran Penduga Parameter beta1 (x1)")
CoordK@data$b2 <- gwr.model$SDF$x2
spplot(CoordK, zcol="b2", main="Peta Sebaran Penduga Parameter beta2 (x2)")
CoordK@data$b3 <- gwr.model$SDF$x3
spplot(CoordK, zcol="b3", main="Peta Sebaran Penduga Parameter beta3 (x3)")
7.2 Geographically Weighted Regression (GWR) Gaussian kernel
Gaussian Kernel:
\[ w_{ij} = \exp \left( -\frac{d_{ij}^2}{2h^2} \right) \]
# Select bandwidth using Gaussian kernel
gwr_gaussian_bw <- gwr.sel(y ~ x1 + x2 + x3, data=data, coords=longlat, gweight=gwr.Gauss)## Bandwidth: 0.894596 CV score: 238105666
## Bandwidth: 1.446042 CV score: 239639540
## Bandwidth: 0.5537839 CV score: 250479142
## Bandwidth: 1.138571 CV score: 238577856
## Bandwidth: 0.7644173 CV score: 239034035
## Bandwidth: 0.9766506 CV score: 238123268
## Bandwidth: 0.9255316 CV score: 238081533
## Bandwidth: 0.9301108 CV score: 238081595
## Bandwidth: 0.9275184 CV score: 238081457
## Bandwidth: 0.9274777 CV score: 238081457
## Bandwidth: 0.927437 CV score: 238081457
## Bandwidth: 0.9274777 CV score: 238081457# Fit GWR model using Gaussian kernel
gwr_gaussian_model <- gwr(y ~ x1 + x2 + x3, data=data, coords=longlat, gweight=gwr.Gauss, bandwidth=gwr_gaussian_bw, hatmatrix=TRUE)
# Display model summary
gwr_gaussian_model## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = data, coords = longlat,
## bandwidth = gwr_gaussian_bw, gweight = gwr.Gauss, hatmatrix = TRUE)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 0.9274777
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -2.5275e+01 2.7974e+02 9.8290e+02 1.9986e+03 2.8266e+03
## x1 1.1910e-01 2.0019e-01 3.0969e-01 4.0699e-01 4.4885e-01
## x2 -1.3299e-02 -1.0405e-02 -7.6170e-03 -3.6534e-03 -1.8173e-03
## x3 2.1055e-01 2.2358e-01 2.4483e-01 2.5777e-01 2.6522e-01
## Global
## X.Intercept. 1105.5696
## x1 0.2761
## x2 -0.0081
## x3 0.2463
## Number of data points: 27
## Effective number of parameters (residual: 2traceS - traceS'S): 7.087318
## Effective degrees of freedom (residual: 2traceS - traceS'S): 19.91268
## Sigma (residual: 2traceS - traceS'S): 2345.384
## Effective number of parameters (model: traceS): 5.817588
## Effective degrees of freedom (model: traceS): 21.18241
## Sigma (model: traceS): 2274.004
## Sigma (ML): 2014.174
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 506.6448
## AIC (GWR p. 96, eq. 4.22): 493.2703
## Residual sum of squares: 109536204
## Quasi-global R2: 0.86887157.3 Koefisien GWR
coef_gwr<-gwr.model$SDF@data[,c("x1","x2","x3")];coef_gwr## x1 x2 x3
## 1 0.06603811 -0.0146046204 0.2676420
## 2 0.11547311 -0.0124897015 0.2628005
## 3 0.24710507 -0.0092223544 0.2543056
## 4 0.32917551 -0.0067978699 0.2426597
## 5 0.39373242 -0.0042995028 0.2316739
## 6 0.47381138 -0.0009571322 0.2094606
## 7 0.49041569 -0.0001171626 0.1987284
## 8 0.48314951 -0.0003770142 0.1972992
## 9 0.46233481 -0.0014037536 0.2004614
## 10 0.41267767 -0.0038540352 0.2186882
## 11 0.36248291 -0.0058840391 0.2329801
## 12 0.36208385 -0.0061944880 0.2281014
## 13 0.28397313 -0.0089568345 0.2466971
## 14 0.23965895 -0.0101637202 0.2539881
## 15 0.17980595 -0.0123415079 0.2601415
## 16 0.11643539 -0.0140393084 0.2658436
## 17 0.27077465 -0.0088645440 0.2508612
## 18 0.51469496 0.0007886955 0.1897404
## 19 0.08009455 -0.0142458063 0.2668638
## 20 0.16293013 -0.0117180331 0.2611603
## 21 0.31296334 -0.0075228897 0.2443257
## 22 0.46243515 -0.0013956542 0.2003639
## 23 0.08650969 -0.0146409409 0.2677177
## 24 0.05705484 -0.0150949683 0.2688296
## 25 0.29305610 -0.0081952361 0.2475082
## 26 0.46806103 -0.0012239015 0.2090549
## 27 0.50569928 0.0006479104 0.1925244Comparison of OLS and GWR models
8.1 using AIC
aic_values <- data.frame(
"MODEL" = c("GWR", "OLS"),
"AIC" = c(gwr.model$results$AICh, AIC(reg.klasik))
)
aic_values## MODEL AIC
## 1 GWR 492.2630
## 2 OLS 503.84548.2 Plotting the actual vs predicted values for GWR and OLS
CoordK@data$pred_gwr <- gwr.model$SDF$pred
CoordK@data$pred_ols <- predict(reg.klasik, data)
spplot(CoordK, c("pred_gwr", "pred_ols", "y"), names.attr=c("Prediksi GWR", "Prediksi OLS", "Nilai Aktual"), as.table=TRUE)
Asumsi GWR
9.1 Autokorelasi
residual_gwr <- gwr.model$SDF$gwr.e
coords <- cbind(data$longitude, data$latitude)
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")
moran_queen <- moran.test(residual_gwr, queen_listw)
print(moran_queen)##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: queen_listw
##
## Moran I statistic standard deviate = -0.068836, p-value = 0.5274
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04809089 -0.03846154 0.019568809.2 uji normality
shapiro.test(residual_gwr)##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.9606, p-value = 0.38159.3 uji heterogenitas
library(GWmodel)
# Ekstrak residual
gwr_residuals <- gwr.model$SDF$gwr.e
# Uji Breusch-Pagan
bp_test <- bptest(gwr_residuals ~ gwr.model$SDF$x1 + gwr.model$SDF$x2 + gwr.model$SDF$x3)
# Tampilkan hasil
print(bp_test)##
## studentized Breusch-Pagan test
##
## data: gwr_residuals ~ gwr.model$SDF$x1 + gwr.model$SDF$x2 + gwr.model$SDF$x3
## BP = 8.3509, df = 3, p-value = 0.039299.5 membuat plot-spplot
coef_x1 <- gwr.model$SDF$x1
coef_x2 <- gwr.model$SDF$x2
coef_x3 <- gwr.model$SDF$x3Signifikansi secara lokal
##Plot koefisien regresi lokal untuk variabel x1##
spplot(gwr.model$SDF, "x1", main = "Koefisien GWR untuk x1", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x2##
spplot(gwr.model$SDF, "x2", main = "Koefisien GWR untuk x2", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x3##
spplot(gwr.model$SDF, "x3", main = "Koefisien GWR untuk x3", col.regions = terrain.colors(100))
# Menghitung residual untuk OLS dan GWR
residual_ols <- residuals(reg.klasik)
residual_gwr <- gwr.model$SDF$gwr.e
# Menghitung Sum of Squares (SS) untuk OLS dan GWR
SS_ols <- sum(residual_ols^2) # Sum of Squares untuk model OLS
SS_gwr <- sum(residual_gwr^2) # Sum of Squares untuk model GWR
# Degrees of Freedom (df) untuk OLS dan GWR
df_ols <- df.residual(reg.klasik) # Derajat bebas model OLS
df_gwr <- sum(gwr.model$results$edf) # Effective degrees of freedom untuk GWR
# Improvement dari GWR terhadap OLS
SS_improvement <- SS_ols - SS_gwr # Perbaikan Sum of Squares
df_improvement <- df_ols - df_gwr # Perbedaan derajat bebas
# Menghitung Mean Square (MS)
MS_ols <- SS_ols / df_ols # Mean Square OLS
MS_improvement <- SS_improvement / df_improvement # Mean Square improvement
MS_gwr <- SS_gwr / df_gwr # Mean Square GWR residual
# Menghitung nilai F-statistik
F_value <- MS_improvement / MS_gwr # F-statistik untuk membandingkan OLS dan GWR
# Menghitung p-value dari F-statistik
p_value <- pf(F_value, df_improvement, df_gwr, lower.tail = FALSE)
# Menampilkan tabel hasil perbandingan OLS dan GWR
tabel_hasil <- data.frame(
"Keragaman Residual" = c("OLS Residual", "GWR Improvement", "GWR Residual"),
"Df" = c(df_ols, df_improvement, df_gwr),
"Sum Sq" = c(SS_ols, SS_improvement, SS_gwr),
"Mean Sq" = c(MS_ols, MS_improvement, MS_gwr),
"F" = c(NA, NA, F_value), # Nilai F hanya muncul untuk GWR Improvement
"p-value" = c(NA, NA, p_value) # Menambahkan kolom p-value
)
# Menampilkan tabel
print(tabel_hasil)## Keragaman.Residual Df Sum.Sq Mean.Sq F p.value
## 1 OLS Residual 23.000000 138797921 6034692 NA NA
## 2 GWR Improvement 3.601862 35127957 9752721 NA NA
## 3 GWR Residual 19.398138 103669964 5344326 1.824874 0.1691947##membuat plot-spplot##
coef_x1 <- gwr.model$SDF$x1
coef_x2 <- gwr.model$SDF$x2
coef_x3 <- gwr.model$SDF$x3
##Membuat plot -gglot2##
gwr_df <- as.data.frame(gwr.model$SDF)# Plot koefisien lokal untuk x1
ggplot(gwr_df, aes(x = data$longitude, y = data$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 = data$longitude, y = data$latitude)) +
geom_point(aes(color = x2), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x2", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x3
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x3), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x3", color = "Koefisien") +
theme_minimal()
# Ambil hasil koefisien dan standar error
hasil_gwr <- as.data.frame(gwr.model$SDF)
# Contoh kolom koefisien dan p-value yang dihasilkan
# Kolom biasanya bernama seperti "x1", "x1_se", "x2", "x2_se", dll.
# Periksa struktur hasil GWR
str(hasil_gwr)## 'data.frame': 27 obs. of 20 variables:
## $ sum.w : num 13 11.8 16.2 19 17 ...
## $ X.Intercept. : num 3589 3043 1570 810 276 ...
## $ x1 : num 0.066 0.115 0.247 0.329 0.394 ...
## $ x2 : num -0.0146 -0.01249 -0.00922 -0.0068 -0.0043 ...
## $ x3 : num 0.268 0.263 0.254 0.243 0.232 ...
## $ X.Intercept._se : num 1516 1549 1112 1001 1032 ...
## $ x1_se : num 0.135 0.136 0.117 0.114 0.119 ...
## $ x2_se : num 0.00528 0.00543 0.00481 0.00464 0.00497 ...
## $ x3_se : num 0.0335 0.0337 0.0318 0.0313 0.0338 ...
## $ gwr.e : num 1519 2132 2747 -1478 -343 ...
## $ pred : num 25393 8601 6233 13773 8744 ...
## $ pred.se : num 1831 738 1233 774 690 ...
## $ localR2 : num 0.845 0.848 0.853 0.857 0.858 ...
## $ X.Intercept._se_EDF: num 1567 1600 1148 1034 1066 ...
## $ x1_se_EDF : num 0.14 0.141 0.121 0.118 0.123 ...
## $ x2_se_EDF : num 0.00546 0.00561 0.00497 0.00479 0.00513 ...
## $ x3_se_EDF : num 0.0346 0.0348 0.0329 0.0323 0.035 ...
## $ pred.se.1 : num 1892 762 1274 800 713 ...
## $ coord.x : num 107 107 107 108 108 ...
## $ coord.y : num -6.56 -7.07 -7.13 -7.1 -7.36 ...T-test for significance of coefficients in GWR
hasil_gwr <- hasil_gwr %>%
mutate(
p_value_x1 = 2 * (1 - pnorm(abs(coef_x1 / x1_se))),
p_value_x2 = 2 * (1 - pnorm(abs(coef_x2 / x2_se))),
p_value_x3 = 2 * (1 - pnorm(abs(coef_x3 / x3_se)))
) %>%
# Tentukan variabel mana yang signifikan (p < 0.05)
mutate(
signifikan_x1 = ifelse(p_value_x1 < 0.05, "x1",NA),
signifikan_x2 = ifelse(p_value_x2 < 0.05, "x2",NA),
signifikan_x3 = ifelse(p_value_x3 < 0.05, "x3",NA)
)# Tambahkan kolom Kabupaten.kota ke hasil_gwr sebelum melakukan analisis
hasil_gwr <- hasil_gwr %>%
mutate(Kabupaten.Kota = data$Kabupaten.Kota) # Pastikan 'data' adalah nama dataframe awal yang berisi Kabupaten.kota# Menggabungkan variabel signifikan menjadi satu kolom
hasil_gwr <- hasil_gwr %>%
mutate(
variabel_signifikan = paste(
signifikan_x1, signifikan_x2, signifikan_x3, sep = ", "
)
) %>%
# Hapus NA dalam hasil variabel_signifikan
mutate(variabel_signifikan = gsub("NA, |, NA", "", variabel_signifikan)) %>%
mutate(variabel_signifikan = ifelse(variabel_signifikan == "NA", "Tidak ada", variabel_signifikan))
# Membuat tabel akhir
tabel_hasil <- hasil_gwr %>%
dplyr::select(Kabupaten.Kota, variabel_signifikan)# Tampilkan tabel hasil
print(tabel_hasil)## Kabupaten.Kota variabel_signifikan
## 1 BOGOR x2, x3
## 2 SUKABUMI x2, x3
## 3 CIANJUR x1, x3
## 4 Bandung x1, x3
## 5 GARUT x1, x3
## 6 TASIKMALAYA x1, x3
## 7 CIAMIS x1, x3
## 8 KUNINGAN x1, x3
## 9 CIREBON x1, x3
## 10 MAJALENGKA x1, x3
## 11 SUMEDANG x1, x3
## 12 INDRAMAYU x1, x3
## 13 SUBANG x1, x3
## 14 Purwakarta x1, x2, x3
## 15 KARAWANG x2, x3
## 16 BEKASI x2, x3
## 17 BANDUNG BARAT x1, x3
## 18 PANGANDARAN x1, x3
## 19 KOTA BOGOR x2, x3
## 20 Kota Sukabumi x2, x3
## 21 KOTA BANDUNG x1, x3
## 22 KOTA CIREBON x1, x3
## 23 KOTA BEKASI x2, x3
## 24 KOTA DEPOK x2, x3
## 25 KOTA CIMAHI x1, x3
## 26 KOTA TASIKMALAYA x1, x3
## 27 KOTA BANJAR x1, x3view(tabel_hasil)