This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
## 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
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
## Loading required package: sp
## corrplot 0.95 loaded
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
##
## Recode
## Loading required package: Matrix
##
## 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.ZeroPolicyOption
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── 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 conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
##
## Loading required package: robustbase
##
## Loading required package: Rcpp
##
## Welcome to GWmodel version 2.4-2.
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
# Set working directory
setwd("C:/Users/suhaila/Downloads")
# Load data
data <- read.csv("spasial - Sheet1 (6).csv", dec=",")
jabar <- read_sf("jabar_modified.shp")
#####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 geometries
ggplot() +
geom_sf(data = jabar_merged, aes(fill = Balita.Gizi.Buruk), color = "white") +
theme_bw() +
scale_fill_gradient(low = "pink", high = "brown") +
# Menambahkan nama kabupaten
geom_sf_text(data = jabar_centroid, aes(label = Kabupaten.kota), size = 3, color = "white") +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank() # Remove y-axis labels
) +
labs(title = "", fill = "Gizi Buruk")
summary(data$Balita.Gizi.Buruk)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 55.0 125.0 242.7 380.5 846.0
#QUEEN
# Menggunakan koordinat dari data untuk menghitung matriks bobot
coords <- cbind(data$longitude, data$latitude)
# 1. Matriks Bobot KNN (contoh 5 dan 6 tetangga terdekat)
k_neighbors1 <- 5
knn_weights1 <- knearneigh(coords, k = k_neighbors1)
knn_nb1 <- knn2nb(knn_weights1)
knn_listw1 <- nb2listw(knn_nb1, style = "W") # Menggunakan bobot standar
k_neighbors2 <- 6
knn_weights2 <- knearneigh(coords, k = k_neighbors2)
knn_nb2 <- knn2nb(knn_weights2)
knn_listw2 <- nb2listw(knn_nb2, style = "W") # Menggunakan bobot standar
# 2. Matriks Bobot Queen
queen_nb <- poly2nb(jabar_sf) # Pastikan jabar_sf adalah data geospasial
queen_listw <- nb2listw(queen_nb, style = "W") # Menggunakan bobot standar
# 3. Matriks Bobot Rook
rook_nb <- poly2nb(jabar_sf, queen = FALSE) # Rook menggunakan queen = FALSE
rook_listw <- nb2listw(rook_nb, style = "W") # Menggunakan bobot standar
# Load library
library(spdep)
# Misalkan 'variabel_target' adalah nama kolom dari data yang ingin diuji
variabel_target <- data$Balita.Gizi.Buruk # Ganti 'nama_kolom_target' sesuai dengan nama kolom Anda
# Hitung nilai Moran's I untuk KNN
knn_moran1 <- moran.test(variabel_target, knn_listw1)
print(knn_moran1)
##
## Moran I test under randomisation
##
## data: variabel_target
## weights: knn_listw1
##
## Moran I statistic standard deviate = -0.86937, p-value = 0.8077
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.123347931 -0.038461538 0.009533903
knn_moran2 <- moran.test(variabel_target, knn_listw2)
print(knn_moran2)
##
## Moran I test under randomisation
##
## data: variabel_target
## weights: knn_listw2
##
## Moran I statistic standard deviate = -1.2181, p-value = 0.8884
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.143005613 -0.038461538 0.007366245
# Hitung nilai Moran's I untuk Queen
queen_moran <- moran.test(variabel_target, queen_listw)
print(queen_moran)
##
## Moran I test under randomisation
##
## data: variabel_target
## weights: queen_listw
##
## Moran I statistic standard deviate = -2.1286, p-value = 0.9834
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.33447439 -0.03846154 0.01933863
# Hitung nilai Moran's I untuk Rook
rook_moran <- moran.test(variabel_target, rook_listw)
print(rook_moran)
##
## Moran I test under randomisation
##
## data: variabel_target
## weights: rook_listw
##
## Moran I statistic standard deviate = -2.1286, p-value = 0.9834
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.33447439 -0.03846154 0.01933863
# Bandingkan hasil Moran's I
moran_results <- data.frame(
Method = c("KNN 5","KNN 6", "Queen", "Rook"),
Moran_I = c(knn_moran1$estimate[1],knn_moran2$estimate[1], queen_moran$estimate[1], rook_moran$estimate[1]),
Expectation = c(knn_moran1$estimate[2],knn_moran2$estimate[2], queen_moran$estimate[2], rook_moran$estimate[2]),
Variance = c(knn_moran1$estimate[3],knn_moran2$estimate[3], queen_moran$estimate[3], rook_moran$estimate[3]),
P_Value = c(knn_moran1$p.value, knn_moran2$p.value, queen_moran$p.value, rook_moran$p.value)
)
print(moran_results)
## Method Moran_I Expectation Variance P_Value
## 1 KNN 5 -0.1233479 -0.03846154 0.009533903 0.8076765
## 2 KNN 6 -0.1430056 -0.03846154 0.007366245 0.8884035
## 3 Queen -0.3344744 -0.03846154 0.019338634 0.9833570
## 4 Rook -0.3344744 -0.03846154 0.019338634 0.9833570
# Daftar variabel untuk dipetakan
vars <- c("PDRB.per.kapita", "tingkat.pengangguran.terbuka",
"Indeks.keparahan.kemiskinan", "kepadatan.penduduk",
"jumlah.dokter", "cakupan.imunisasi.dasar",
"keluarga.dengan.akses.fasilitas.sanitasi",
"balita.berat.badan.kurang")
# Loop untuk membuat peta setiap variabel
for (var in vars) {
# Plot untuk setiap variabel
p <- ggplot() +
geom_sf(data = jabar_merged, aes_string(fill = var), color = "white") +
scale_fill_gradient(low = "lightblue", high = "darkblue", na.value = "grey90") +
theme_bw() +
geom_sf_text(data = jabar_centroid, aes(label = Kabupaten.kota),
size = 3, color = "white", fontface = "bold") +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right",
axis.text.x = element_blank(),
axis.text.y = element_blank()
) +
labs(title = paste("Peta Distribusi", var), fill = var)
# Menampilkan plot
print(p)
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Deskriptif Data
summary(data)
## Kabupaten.kota Balita.Gizi.Buruk PDRB.per.kapita
## Length:27 Min. : 0.0 Min. :14205
## Class :character 1st Qu.: 55.0 1st Qu.:18882
## Mode :character Median :125.0 Median :23973
## Mean :242.7 Mean :32041
## 3rd Qu.:380.5 3rd Qu.:35051
## Max. :846.0 Max. :88554
## tingkat.pengangguran.terbuka Indeks.keparahan.kemiskinan kepadatan.penduduk
## Min. : 1.520 Min. :0.0300 Min. : 382.0
## 1st Qu.: 6.535 1st Qu.:0.2000 1st Qu.: 824.5
## Median : 7.650 Median :0.2600 Median : 1449.0
## Mean : 7.186 Mean :0.2807 Mean : 3873.3
## 3rd Qu.: 8.500 3rd Qu.:0.3450 3rd Qu.: 5749.0
## Max. :10.520 Max. :0.6200 Max. :15047.0
## jumlah.puskesmas jumlah.dokter cakupan.imunisasi.dasar
## Min. : 10.00 Min. : 45.0 Min. : 82.67
## 1st Qu.: 23.50 1st Qu.: 76.0 1st Qu.: 95.51
## Median : 38.00 Median : 96.0 Median : 97.14
## Mean : 40.96 Mean :114.1 Mean : 97.52
## 3rd Qu.: 50.50 3rd Qu.:135.5 3rd Qu.: 99.67
## Max. :101.00 Max. :237.0 Max. :123.83
## keluarga.dengan.akses.fasilitas.sanitasi
## Min. : 54951
## 1st Qu.: 221323
## Median : 468556
## Mean : 481529
## 3rd Qu.: 625727
## Max. :1119445
## jumlah.penduduk.dengan.akses.berkelanjutan.air.minum.layak
## Min. : 182879
## 1st Qu.: 917626
## Median :1164650
## Mean :1422445
## 3rd Qu.:2026029
## Max. :2966480
## balita.berat.badan.kurang longitude latitude
## Min. : 638 Min. :-7701397 Min. : 1081971
## 1st Qu.: 2482 1st Qu.:-6994542 1st Qu.:106878230
## Median : 4482 Median :-6836531 Median :107512302
## Mean : 5845 Mean :-6605407 Mean : 92565328
## 3rd Qu.: 6771 3rd Qu.:-6515690 3rd Qu.:108016295
## Max. :15799 Max. : -687121 Max. :108557818
## id
## Min. : 1.0
## 1st Qu.: 7.5
## Median :14.0
## Mean :14.0
## 3rd Qu.:20.5
## Max. :27.0
#OLS#
y <- as.numeric(data$Balita.Gizi.Buruk)
x1 <- as.numeric(data$PDRB.per.kapita)
x2 <- data$tingkat.pengangguran.terbuka
x3 <- data$Indeks.keparahan.kemiskinan
x4 <- as.numeric(data$kepadatan.penduduk)
x5 <- data$jumlah.puskesmas
x6 <- data$jumlah.dokter
x7 <- data$cakupan.imunisasi.dasar
x8 <- data$keluarga.dengan.akses.fasilitas.sanitasi
x9 <- data$jumlah.penduduk.dengan.akses.berkelanjutan
x10 <- data$balita.berat.badan.kurang
model_ols <- lm(y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 +x10, data = data)
summary(model_ols)
##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x10, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -211.60 -65.80 -32.12 71.38 215.83
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.022e+02 5.013e+02 -1.800 0.08871 .
## x1 2.510e-03 1.633e-03 1.537 0.14173
## x2 -3.486e+01 1.693e+01 -2.059 0.05424 .
## x3 7.271e+02 2.447e+02 2.971 0.00819 **
## x4 1.170e-02 9.355e-03 1.251 0.22693
## x6 1.229e+00 7.466e-01 1.646 0.11709
## x7 6.988e+00 5.213e+00 1.340 0.19677
## x8 2.958e-04 1.639e-04 1.804 0.08794 .
## x10 1.735e-02 1.199e-02 1.446 0.16523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 136 on 18 degrees of freedom
## Multiple R-squared: 0.7926, Adjusted R-squared: 0.7004
## F-statistic: 8.597 on 8 and 18 DF, p-value: 8.426e-05
model_ols1 <- lm(y ~ x3, data = data)
summary(model_ols1)
##
## Call:
## lm(formula = y ~ x3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -285.7 -173.2 -106.7 144.3 571.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 59.37 103.37 0.574 0.5709
## x3 653.19 330.88 1.974 0.0595 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 235.6 on 25 degrees of freedom
## Multiple R-squared: 0.1349, Adjusted R-squared: 0.1003
## F-statistic: 3.897 on 1 and 25 DF, p-value: 0.05952
##Heterogenitas##
bptest(model_ols)
##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 16.221, df = 8, p-value = 0.03932
##Multikolinearity##
vif(model_ols)
## x1 x2 x3 x4 x6 x7 x8 x10
## 1.693171 1.649316 1.643057 2.633562 2.241785 2.098027 3.527768 3.766920
##Normalitas##
residual_ols <- residuals(model_ols)
shapiro.test(residual_ols)
##
## Shapiro-Wilk normality test
##
## data: residual_ols
## W = 0.95746, p-value = 0.3229
##Autokorelasi##
coords <- cbind(data$longitude, data$latitude)
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")
moran.test(residual_ols, queen_listw)
##
## Moran I test under randomisation
##
## data: residual_ols
## weights: queen_listw
##
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.03355745 -0.03846154 0.01968442
# Hitung nilai Moran's I untuk KNN
nb_5 <- knn2nb(knearneigh(coords, k=5))
listw_5 <- nb2listw(nb_5, style="W")
nb_6 <- knn2nb(knearneigh(coords, k=6))
listw_6 <- nb2listw(nb_6, style="W")
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")
rook_nb <- poly2nb(jabar_sf, queen = FALSE)
rook_listw <- nb2listw(rook_nb, style = "W")
knn_moran_5 <- moran.test(residual_ols, listw_5)
print(knn_moran_5)
##
## Moran I test under randomisation
##
## data: residual_ols
## weights: listw_5
##
## Moran I statistic standard deviate = 1.251, p-value = 0.1055
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.084773021 -0.038461538 0.009703921
knn_moran_6 <- moran.test(residual_ols, listw_6)
print(knn_moran_6)
##
## Moran I test under randomisation
##
## data: residual_ols
## weights: listw_6
##
## Moran I statistic standard deviate = 1.8443, p-value = 0.03257
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.121219928 -0.038461538 0.007496251
knn_moran_7 <- moran.test(residual_ols, queen_listw)
print(knn_moran_7)
##
## Moran I test under randomisation
##
## data: residual_ols
## weights: queen_listw
##
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.03355745 -0.03846154 0.01968442
knn_moran_8 <- moran.test(residual_ols, rook_listw)
print(knn_moran_8)
##
## Moran I test under randomisation
##
## data: residual_ols
## weights: rook_listw
##
## Moran I statistic standard deviate = 0.51332, p-value = 0.3039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.03355745 -0.03846154 0.01968442
######GWR#####
y <- as.numeric(data$Balita.Gizi.Buruk)
x1 <- as.numeric(data$PDRB.per.kapita)
x2 <- data$tingkat.pengangguran.terbuka
x3 <- data$Indeks.keparahan.kemiskinan
x4 <- as.numeric(data$kepadatan.penduduk)
x5 <- data$jumlah.puskesmas
x6 <- data$jumlah.dokter
x7 <- data$cakupan.imunisasi.dasar
x8 <- data$keluarga.dengan.akses.fasilitas.sanitasi
x9 <- data$jumlah.penduduk.dengan.akses.berkelanjutan
x10 <- data$balita.berat.badan.kurang
#optimal bandwith
bandwidth_optimal <- gwr.sel(
y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 +x10,
data = data,
coords = cbind(data$longitude, data$latitude),
adapt = TRUE
)
## Adaptive q: 0.381966 CV score: 801523.8
## Adaptive q: 0.618034 CV score: 736303.4
## Adaptive q: 0.763932 CV score: 738355.4
## Adaptive q: 0.6817313 CV score: 736847.3
## Adaptive q: 0.586383 CV score: 740832.3
## Adaptive q: 0.6180747 CV score: 736297.9
## Adaptive q: 0.6423894 CV score: 734862.3
## Adaptive q: 0.639664 CV score: 734853.9
## Adaptive q: 0.6404948 CV score: 734855.9
## Adaptive q: 0.6314176 CV score: 734867.8
## Adaptive q: 0.6371811 CV score: 734851.6
## Adaptive q: 0.6373114 CV score: 734851.6
## Adaptive q: 0.6373521 CV score: 734851.6
## Adaptive q: 0.6373928 CV score: 734851.6
## Adaptive q: 0.6373521 CV score: 734851.6
#running GWR
model_gwr <- gwr(
y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 +x10,
data = data,
coords = cbind(data$longitude, data$latitude),
adapt = bandwidth_optimal,
hatmatrix = TRUE,
se.fit = TRUE
)
print(model_gwr)
## Call:
## gwr(formula = y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x10, data = data,
## coords = cbind(data$longitude, data$latitude), adapt = bandwidth_optimal,
## hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.6373521 (about 17 of 27 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. -1.2595e+03 -1.2179e+03 -1.1656e+03 -1.0903e+03 -8.1918e+02
## x1 2.4099e-03 2.5537e-03 2.6297e-03 2.7340e-03 2.8809e-03
## x2 -3.9520e+01 -3.8227e+01 -3.5289e+01 -3.2945e+01 -3.2371e+01
## x3 6.5091e+02 8.4012e+02 8.7225e+02 1.0354e+03 1.0736e+03
## x4 1.0290e-02 1.3426e-02 1.6364e-02 2.0355e-02 2.1841e-02
## x6 1.0784e+00 1.1317e+00 1.3719e+00 1.5923e+00 1.6265e+00
## x7 6.6279e+00 7.3115e+00 8.2238e+00 9.0969e+00 9.6672e+00
## x8 3.0045e-04 3.0598e-04 3.1884e-04 3.4699e-04 3.6578e-04
## x10 1.4703e-02 1.6071e-02 1.7021e-02 1.7348e-02 1.8247e-02
## Global
## X.Intercept. -902.1530
## x1 0.0025
## x2 -34.8597
## x3 727.0945
## x4 0.0117
## x6 1.2290
## x7 6.9877
## x8 0.0003
## x10 0.0173
## Number of data points: 27
## Effective number of parameters (residual: 2traceS - traceS'S): 13.68578
## Effective degrees of freedom (residual: 2traceS - traceS'S): 13.31422
## Sigma (residual: 2traceS - traceS'S): 134.4216
## Effective number of parameters (model: traceS): 12.38756
## Effective degrees of freedom (model: traceS): 14.61244
## Sigma (model: traceS): 128.3115
## Sigma (ML): 94.39412
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 379.5052
## AIC (GWR p. 96, eq. 4.22): 334.5741
## Residual sum of squares: 240576.8
## Quasi-global R2: 0.8500062
#koefisien pada setiap lokasi
coef_gwr <- model_gwr$SDF@data[, c("x1",'x2','x3','x4','x6','x7','x8','x10')]
# Menampilkan ringkasan model
summary(model_gwr)
## Length Class Mode
## SDF 27 SpatialPointsDataFrame S4
## lhat 729 -none- numeric
## lm 11 -none- list
## results 14 -none- list
## bandwidth 27 -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
# Menampilkan ringkasan GWR Global (kuartil, median, dll.)
summary_gwr <- apply(coef_gwr, 2, summary)
print(summary_gwr)
## x1 x2 x3 x4 x6 x7
## Min. 0.002409909 -39.51971 650.9142 0.01029028 1.078432 6.627918
## 1st Qu. 0.002553685 -38.22654 840.1179 0.01342622 1.131732 7.311494
## Median 0.002629708 -35.28912 872.2458 0.01636355 1.371947 8.223832
## Mean 0.002628406 -35.81804 896.4846 0.01628267 1.346413 8.191521
## 3rd Qu. 0.002733953 -32.94530 1035.4412 0.02035543 1.592290 9.096944
## Max. 0.002880885 -32.37063 1073.6475 0.02184072 1.626486 9.667184
## x8 x10
## Min. 0.0003004484 0.01470271
## 1st Qu. 0.0003059780 0.01607122
## Median 0.0003188450 0.01702093
## Mean 0.0003257273 0.01680377
## 3rd Qu. 0.0003469948 0.01734795
## Max. 0.0003657821 0.01824698
# Asumsi Heterogenitas Spasial
# Memasang paket lmtest
library(lmtest)
# Melakukan uji Breusch-Pagan pada model GWR
bp_test <- bptest(model_gwr$SDF$gwr.e ~ model_gwr$SDF$pred +
I(model_gwr$SDF$pred^2))
# Menampilkan hasil uji
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: model_gwr$SDF$gwr.e ~ model_gwr$SDF$pred + I(model_gwr$SDF$pred^2)
## BP = 3.9297, df = 2, p-value = 0.1402
#Asumsi Multikoleniaritas
# Mengakses Condition Number dari model GWR
condition_number <- model_gwr$SDF@data$CN
print(condition_number)
## NULL
# Menyaring lokasi dengan potensi multikolinearitas tinggi
if (any(condition_number >= 30)) {
print("Terdapat indikasi multikolinearitas di beberapa lokasi.")
} else {
print("Tidak ada indikasi multikolinearitas.")
}
## [1] "Tidak ada indikasi multikolinearitas."
##uji normality##
residual_gwr <- model_gwr$SDF$gwr.e
shapiro.test(residual_gwr)
##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.94385, p-value = 0.1516
##Autokorelasi##
residual_gwr <- model_gwr$SDF$gwr.e
coords <- cbind(data$longitude, data$latitude)
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")
moran.test(residual_gwr, queen_listw)
##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: queen_listw
##
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12384698 -0.03846154 0.01932658
# Hitung nilai Moran's I untuk KNN
nb5 <- knn2nb(knearneigh(coords, k=5))
listw5 <- nb2listw(nb5, style="W")
nb6 <- knn2nb(knearneigh(coords, k=6))
listw6 <- nb2listw(nb6, style="W")
queen_nb <- poly2nb(jabar_sf)
queen_listw <- nb2listw(queen_nb, style = "W")
rook_nb <- poly2nb(jabar_sf, queen = FALSE)
rook_listw <- nb2listw(rook_nb, style = "W")
knn_moran5 <- moran.test(residual_gwr, listw5)
print(knn_moran5)
##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: listw5
##
## Moran I statistic standard deviate = 1.6641, p-value = 0.04805
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.123973735 -0.038461538 0.009527975
knn_moran6 <- moran.test(residual_gwr, listw6)
print(knn_moran6)
##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: listw6
##
## Moran I statistic standard deviate = 2.2181, p-value = 0.01328
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.151848124 -0.038461538 0.007361712
knn_moran7 <- moran.test(residual_gwr, queen_listw)
print(knn_moran7)
##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: queen_listw
##
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12384698 -0.03846154 0.01932658
knn_moran8 <- moran.test(residual_gwr, rook_listw)
print(knn_moran8)
##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: rook_listw
##
## Moran I statistic standard deviate = 1.1675, p-value = 0.1215
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.12384698 -0.03846154 0.01932658
# Menghitung residual untuk OLS dan GWR
residual_ols <- residuals(model_ols)
residual_gwr <- model_gwr$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(model_ols) # Derajat bebas model OLS
df_gwr <- sum(model_gwr$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 18.00000 332689.22 18482.73 NA NA
## 2 GWR Improvement 4.68578 92112.47 19657.87 NA NA
## 3 GWR Residual 13.31422 240576.75 18069.16 1.087924 0.4084436
# Comparison of OLS and GWR models using AIC
aic_values <- data.frame(
"MODEL" = c("GWR", "OLS"),
"AIC" = c(model_gwr$results$AICh, AIC(model_ols))
)
aic_values
## MODEL AIC
## 1 GWR 334.5741
## 2 OLS 350.9391
# Menampilkan koefisien dan p-value
coef_table <- as.data.frame(model_gwr$SDF@data)
print(coef_table[, c("x1", "x2", "x3", "x4", "x6", "x7", "x8", "x10")])
## x1 x2 x3 x4 x6 x7 x8
## 1 0.002581543 -32.42267 842.1562 0.01345003 1.596316 8.985352 0.0003114371
## 2 0.002516878 -32.79590 842.4880 0.01400564 1.592530 9.256080 0.0003187514
## 3 0.002541646 -33.46475 834.4039 0.01412889 1.613452 9.483595 0.0003148354
## 4 0.002581912 -37.06315 916.7289 0.01767266 1.435713 9.128609 0.0003347314
## 5 0.002605647 -39.51971 1054.3479 0.02184072 1.130564 7.845026 0.0003617093
## 6 0.002620913 -38.79672 1073.6475 0.02179958 1.100690 7.398352 0.0003657821
## 7 0.002411897 -34.76915 651.4759 0.01031201 1.078434 6.627919 0.0003008658
## 8 0.002753446 -38.40670 1060.3660 0.02089069 1.132900 7.271602 0.0003546347
## 9 0.002791675 -38.26755 1049.1771 0.02043927 1.158007 7.323131 0.0003496963
## 10 0.002780088 -38.79380 1049.6114 0.02079337 1.151643 7.483515 0.0003501047
## 11 0.002761220 -39.50682 1024.1309 0.02070760 1.196582 8.034534 0.0003452199
## 12 0.002880885 -38.18553 1023.1510 0.01948506 1.216365 7.493316 0.0003375968
## 13 0.002845187 -38.59977 949.1596 0.01824999 1.371947 8.662774 0.0003225949
## 14 0.002690965 -35.28912 846.0544 0.01478916 1.604036 9.480896 0.0003059643
## 15 0.002409909 -34.75589 650.9142 0.01029028 1.080224 6.629720 0.0003004557
## 16 0.002637363 -33.09434 831.7497 0.01340241 1.626486 9.191968 0.0003037647
## 17 0.002632095 -37.19573 872.2458 0.01665451 1.533344 9.667184 0.0003197521
## 18 0.002636108 -37.83090 1069.2345 0.02111588 1.119273 7.188613 0.0003653089
## 19 0.002565723 -32.42859 844.1474 0.01357196 1.591087 9.001900 0.0003136732
## 20 0.002516569 -32.79627 842.7081 0.01401161 1.592050 9.253847 0.0003188450
## 21 0.002649549 -37.94518 928.5168 0.01815920 1.409497 9.065279 0.0003321409
## 22 0.002798780 -38.11182 1046.7515 0.02027158 1.164664 7.299856 0.0003487698
## 23 0.002629708 -32.62007 839.6238 0.01328203 1.606199 8.949836 0.0003059917
## 24 0.002590905 -32.37063 840.6120 0.01334146 1.599871 8.966540 0.0003099220
## 25 0.002714460 -32.53098 919.2182 0.01636355 1.492462 8.223832 0.0003007749
## 26 0.002411896 -34.76915 651.4753 0.01031201 1.078432 6.627918 0.0003008659
## 27 0.002409990 -34.75622 650.9873 0.01029099 1.080393 6.629882 0.0003004484
## x10
## 1 0.01698135
## 2 0.01599085
## 3 0.01572963
## 4 0.01485624
## 5 0.01604837
## 6 0.01683812
## 7 0.01718762
## 8 0.01779161
## 9 0.01794844
## 10 0.01758856
## 11 0.01664293
## 12 0.01824698
## 13 0.01659912
## 14 0.01609408
## 15 0.01718457
## 16 0.01702093
## 17 0.01470271
## 18 0.01732378
## 19 0.01680949
## 20 0.01598902
## 21 0.01511680
## 22 0.01804406
## 23 0.01737211
## 24 0.01711993
## 25 0.01810237
## 26 0.01718762
## 27 0.01718450
##membuat plot-spplot##
coef_x1 <- model_gwr$SDF$x1
coef_x2 <- model_gwr$SDF$x2
coef_x3 <- model_gwr$SDF$x3
coef_x4 <- model_gwr$SDF$x4
coef_x5 <- model_gwr$SDF$x5
coef_x6 <- model_gwr$SDF$x6
coef_x7 <- model_gwr$SDF$x7
coef_x8 <- model_gwr$SDF$x8
coef_x9 <- model_gwr$SDF$x9
coef_x10 <- model_gwr$SDF$x10
##Membuat plot -gglot2##
gwr_df <- as.data.frame(model_gwr$SDF)
gwr_df
## sum.w X.Intercept. x1 x2 x3 x4 x6
## 1 17.37470 -1213.6388 0.002581543 -32.42267 842.1562 0.01345003 1.596316
## 2 17.24262 -1235.3255 0.002516878 -32.79590 842.4880 0.01400564 1.592530
## 3 17.28654 -1251.5650 0.002541646 -33.46475 834.4039 0.01412889 1.613452
## 4 16.57260 -1220.6626 0.002581912 -37.06315 916.7289 0.01767266 1.435713
## 5 16.04412 -1127.8756 0.002605647 -39.51971 1054.3479 0.02184072 1.130564
## 6 16.51345 -1098.2326 0.002620913 -38.79672 1073.6475 0.02179958 1.100690
## 7 17.95773 -819.2625 0.002411897 -34.76915 651.4759 0.01031201 1.078434
## 8 17.40602 -1088.2575 0.002753446 -38.40670 1060.3660 0.02089069 1.132900
## 9 17.51910 -1091.2367 0.002791675 -38.26755 1049.1771 0.02043927 1.158007
## 10 17.54177 -1102.5894 0.002780088 -38.79380 1049.6114 0.02079337 1.151643
## 11 16.89933 -1140.4684 0.002761220 -39.50682 1024.1309 0.02070760 1.196582
## 12 16.68149 -1101.4723 0.002880885 -38.18553 1023.1510 0.01948506 1.216365
## 13 16.07198 -1185.2137 0.002845187 -38.59977 949.1596 0.01824999 1.371947
## 14 16.46443 -1249.2509 0.002690965 -35.28912 846.0544 0.01478916 1.604036
## 15 17.95278 -819.1775 0.002409909 -34.75589 650.9142 0.01029028 1.080224
## 16 16.90314 -1228.8256 0.002637363 -33.09434 831.7497 0.01340241 1.626486
## 17 16.10434 -1259.4827 0.002632095 -37.19573 872.2458 0.01665451 1.533344
## 18 16.25037 -1084.1768 0.002636108 -37.83090 1069.2345 0.02111588 1.119273
## 19 17.45931 -1215.2295 0.002565723 -32.42859 844.1474 0.01357196 1.591087
## 20 17.24052 -1235.1637 0.002516569 -32.79627 842.7081 0.01401161 1.592050
## 21 16.52328 -1214.6486 0.002649549 -37.94518 928.5168 0.01815920 1.409497
## 22 17.44059 -1089.4467 0.002798780 -38.11182 1046.7515 0.02027158 1.164664
## 23 16.90283 -1210.1162 0.002629708 -32.62007 839.6238 0.01328203 1.606199
## 24 17.20623 -1212.0194 0.002590905 -32.37063 840.6120 0.01334146 1.599871
## 25 14.61569 -1165.6284 0.002714460 -32.53098 919.2182 0.01636355 1.492462
## 26 17.95767 -819.2620 0.002411896 -34.76915 651.4753 0.01031201 1.078432
## 27 17.95824 -819.2357 0.002409990 -34.75622 650.9873 0.01029099 1.080393
## x7 x8 x10 X.Intercept._se x1_se x2_se
## 1 8.985352 0.0003114371 0.01698135 517.6940 0.001766681 17.64267
## 2 9.256080 0.0003187514 0.01599085 517.8146 0.001768092 17.59790
## 3 9.483595 0.0003148354 0.01572963 518.2396 0.001771664 17.77064
## 4 9.128609 0.0003347314 0.01485624 514.4076 0.001770659 17.60831
## 5 7.845026 0.0003617093 0.01604837 516.9336 0.001771186 17.42070
## 6 7.398352 0.0003657821 0.01683812 516.4236 0.001766353 17.36570
## 7 6.627919 0.0003008658 0.01718762 476.8166 0.001565843 16.12392
## 8 7.271602 0.0003546347 0.01779161 514.3254 0.001761514 17.36620
## 9 7.323131 0.0003496963 0.01794844 513.8364 0.001759948 17.37226
## 10 7.483515 0.0003501047 0.01758856 514.1727 0.001761698 17.40229
## 11 8.034534 0.0003452199 0.01664293 514.3090 0.001766845 17.53687
## 12 7.493316 0.0003375968 0.01824698 513.6876 0.001758481 17.47963
## 13 8.662774 0.0003225949 0.01659912 514.1146 0.001770960 18.15236
## 14 9.480896 0.0003059643 0.01609408 517.1864 0.001777348 18.43094
## 15 6.629720 0.0003004557 0.01718457 476.7919 0.001566241 16.12237
## 16 9.191968 0.0003037647 0.01702093 519.0886 0.001771919 18.05228
## 17 9.667184 0.0003197521 0.01470271 515.4180 0.001779836 18.11496
## 18 7.188613 0.0003653089 0.01732378 515.3387 0.001762553 17.33894
## 19 9.001900 0.0003136732 0.01680949 517.4553 0.001766254 17.59931
## 20 9.253847 0.0003188450 0.01598902 517.7968 0.001768050 17.59594
## 21 9.065279 0.0003321409 0.01511680 513.9146 0.001772232 17.73193
## 22 7.299856 0.0003487698 0.01804406 513.6803 0.001759406 17.36671
## 23 8.949836 0.0003059917 0.01737211 518.3167 0.001768187 17.81738
## 24 8.966540 0.0003099220 0.01711993 517.9354 0.001766975 17.67191
## 25 8.223832 0.0003007749 0.01810237 502.1377 0.001752437 16.92633
## 26 6.627918 0.0003008659 0.01718762 476.8166 0.001565843 16.12392
## 27 6.629882 0.0003004484 0.01718450 476.7852 0.001566194 16.12203
## x3_se x4_se x6_se x7_se x8_se x10_se gwr.e
## 1 280.1248 0.011785217 1.0900669 5.246787 0.0002234409 0.01150636 175.563911
## 2 281.7249 0.011796945 1.0905776 5.244681 0.0002233612 0.01150831 224.985089
## 3 282.4750 0.011815502 1.0924069 5.256962 0.0002238421 0.01153569 10.531596
## 4 280.4130 0.011832863 1.0885846 5.202321 0.0002246551 0.01145189 -157.928870
## 5 280.1940 0.011990876 1.1097631 5.221019 0.0002278874 0.01148939 -136.410040
## 6 279.2165 0.011928970 1.1058535 5.216955 0.0002260751 0.01146668 183.868125
## 7 234.2089 0.008919521 0.7172558 4.934976 0.0001560973 0.01131969 51.455767
## 8 277.7656 0.011838650 1.0992553 5.194804 0.0002238297 0.01142352 -19.166125
## 9 277.7068 0.011809685 1.0961395 5.190393 0.0002232762 0.01140951 13.920179
## 10 277.8274 0.011841943 1.1007064 5.195290 0.0002239785 0.01142092 -18.162464
## 11 278.0422 0.011876261 1.1041305 5.202932 0.0002253272 0.01144142 -99.103542
## 12 278.8616 0.011774739 1.0931372 5.197286 0.0002228017 0.01140405 32.977222
## 13 279.7379 0.011798742 1.0989727 5.236319 0.0002249652 0.01148187 -77.443513
## 14 281.9936 0.011820438 1.0967432 5.276254 0.0002252315 0.01157563 -47.122907
## 15 234.2509 0.008920980 0.7169945 4.934879 0.0001560783 0.01131972 -3.426360
## 16 281.2947 0.011813487 1.0943764 5.279754 0.0002244971 0.01156964 -24.468625
## 17 282.5514 0.011864470 1.0958203 5.240003 0.0002257659 0.01152628 -41.515225
## 18 278.1688 0.011846946 1.0968161 5.206380 0.0002242518 0.01143639 -57.438057
## 19 280.1769 0.011781869 1.0896008 5.242031 0.0002232933 0.01149838 -30.103397
## 20 281.7158 0.011796694 1.0905404 5.244387 0.0002233540 0.01150779 102.358368
## 21 279.5843 0.011834147 1.0915713 5.204729 0.0002251158 0.01145524 -9.654120
## 22 277.6673 0.011798266 1.0945212 5.188585 0.0002230461 0.01140520 60.359725
## 23 280.2147 0.011796414 1.0916385 5.261465 0.0002239315 0.01153236 -56.977484
## 24 280.1342 0.011788656 1.0905096 5.250839 0.0002235630 0.01151303 -1.071166
## 25 268.5842 0.011665069 1.0540869 5.048940 0.0002135093 0.01133336 30.839991
## 26 234.2090 0.008919522 0.7172560 4.934976 0.0001560973 0.01131969 -86.364383
## 27 234.2456 0.008920811 0.7169707 4.934855 0.0001560762 0.01131972 163.070077
## pred pred.se localR2 X.Intercept._se_EDF x1_se_EDF x2_se_EDF
## 1 592.436089 86.56403 0.8470221 542.3462 0.001850808 18.48280
## 2 382.014911 57.04700 0.8392848 542.4725 0.001852288 18.43590
## 3 103.468404 53.38254 0.8384805 542.9178 0.001856030 18.61686
## 4 641.928870 75.16825 0.8362265 538.9033 0.001854976 18.44680
## 5 281.410040 81.03584 0.8488398 541.5496 0.001855529 18.25026
## 6 346.131875 72.46992 0.8575739 541.0153 0.001850466 18.19264
## 7 122.544233 63.01958 0.8390852 499.5222 0.001640407 16.89173
## 8 67.166125 70.57817 0.8698323 538.8172 0.001845396 18.19316
## 9 832.079821 121.01059 0.8718479 538.3049 0.001843755 18.19951
## 10 124.162464 59.67010 0.8687081 538.6571 0.001845588 18.23097
## 11 157.103542 44.93975 0.8591907 538.8000 0.001850980 18.37196
## 12 585.022778 88.90498 0.8758457 538.1490 0.001842219 18.31200
## 13 363.443513 88.69399 0.8596381 538.5963 0.001855291 19.01676
## 14 112.122907 67.87514 0.8467898 541.8144 0.001861984 19.30861
## 15 241.426360 83.02261 0.8390021 499.4963 0.001640825 16.89011
## 16 349.468625 112.19045 0.8490439 543.8071 0.001856296 18.91191
## 17 166.515225 43.93266 0.8361481 539.9618 0.001864591 18.97758
## 18 67.438057 105.07852 0.8623499 539.8788 0.001846484 18.16460
## 19 30.103397 78.70999 0.8456103 542.0961 0.001850362 18.43738
## 20 -57.358368 65.91062 0.8392705 542.4538 0.001852243 18.43385
## 21 445.654120 105.50254 0.8407535 538.3869 0.001856625 18.57631
## 22 -8.359725 86.77570 0.8726698 538.1414 0.001843188 18.19370
## 23 312.977484 94.61529 0.8507495 542.9985 0.001852386 18.66583
## 24 59.071166 106.11628 0.8480191 542.5991 0.001851117 18.51344
## 25 -12.839991 81.32986 0.8723448 526.0492 0.001835886 17.73235
## 26 204.364383 94.33970 0.8390851 499.5223 0.001640407 16.89173
## 27 -139.070077 67.63229 0.8390076 499.4893 0.001640775 16.88975
## x3_se_EDF x4_se_EDF x6_se_EDF x7_se_EDF x8_se_EDF x10_se_EDF pred.se.1
## 1 293.4641 0.012346420 1.1419750 5.496635 0.0002340809 0.01205428 90.68614
## 2 295.1404 0.012358706 1.1425100 5.494428 0.0002339974 0.01205633 59.76353
## 3 295.9262 0.012378147 1.1444265 5.507294 0.0002345013 0.01208501 55.92458
## 4 293.7661 0.012396335 1.1404221 5.450051 0.0002353530 0.01199722 78.74770
## 5 293.5366 0.012561872 1.1626091 5.469640 0.0002387392 0.01203650 84.89471
## 6 292.5125 0.012497018 1.1585134 5.465382 0.0002368406 0.01201271 75.92088
## 7 245.3618 0.009344261 0.7514109 5.169975 0.0001635305 0.01185872 66.02052
## 8 290.9926 0.012402397 1.1516009 5.442176 0.0002344883 0.01196750 73.93905
## 9 290.9309 0.012372053 1.1483368 5.437555 0.0002339085 0.01195283 126.77302
## 10 291.0573 0.012405846 1.1531212 5.442686 0.0002346442 0.01196478 62.51155
## 11 291.2823 0.012441799 1.1567083 5.450691 0.0002360571 0.01198625 47.07975
## 12 292.1407 0.012335443 1.1451915 5.444777 0.0002334114 0.01194710 93.13857
## 13 293.0587 0.012360589 1.1513049 5.485668 0.0002356779 0.01202863 92.91753
## 14 295.4219 0.012383318 1.1489692 5.527505 0.0002359568 0.01212685 71.10730
## 15 245.4057 0.009345790 0.7511372 5.169874 0.0001635106 0.01185876 86.97608
## 16 294.6897 0.012376036 1.1464897 5.531171 0.0002351874 0.01212057 117.53287
## 17 296.0063 0.012429446 1.1480024 5.489527 0.0002365166 0.01207515 46.02470
## 18 291.4150 0.012411088 1.1490456 5.454304 0.0002349305 0.01198098 110.08227
## 19 293.5187 0.012342913 1.1414867 5.491652 0.0002339264 0.01204592 82.45809
## 20 295.1309 0.012358443 1.1424710 5.494121 0.0002339899 0.01205578 69.04923
## 21 292.8979 0.012397680 1.1435510 5.452574 0.0002358356 0.01200073 110.52648
## 22 290.8896 0.012360091 1.1466415 5.435662 0.0002336673 0.01194830 90.90789
## 23 293.5583 0.012358150 1.1436215 5.512012 0.0002345949 0.01208152 99.12079
## 24 293.4740 0.012350023 1.1424388 5.500880 0.0002342089 0.01206127 111.16946
## 25 281.3740 0.012220550 1.1042817 5.289366 0.0002236765 0.01187305 85.20272
## 26 245.3618 0.009344262 0.7514111 5.169976 0.0001635305 0.01185872 98.83208
## 27 245.4002 0.009345613 0.7511123 5.169849 0.0001635084 0.01185875 70.85289
## coord.x coord.y
## 1 -6479679 106824965
## 2 -6915727 106932576
## 3 -6822558 107139542
## 4 -7012851 107528627
## 5 -7202988 107885592
## 6 -7361212 108112488
## 7 -7325788 1083514
## 8 -6976233 108482982
## 9 -6764507 108478858
## 10 -6836531 108227876
## 11 -6861126 107920102
## 12 -6327234 108322903
## 13 -6571549 107762495
## 14 -6551701 107446541
## 15 -6301721 10730529
## 16 -6364614 107172509
## 17 -6840651 107512302
## 18 -7701397 108495155
## 19 -6594946 106794913
## 20 -6918366 106931496
## 21 -6910786 107609757
## 22 -6707076 108557818
## 23 -6236221 106994293
## 24 -6394473 106822692
## 25 -687121 107555486
## 26 -7316436 1081971
## 27 -7362487 10855887
jabar_merged1 <- jabar_merged %>%
bind_cols(gwr_df %>%
select(X.Intercept., x1, x2, x3, x4, x6, x7, x8, x10))
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = X.Intercept.), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk intercept", fill = "Koefisien Beta0") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x1), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x1", fill = "Koefisien Beta1") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x2), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x2", fill = "Koefisien Beta2") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x3), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x3", fill = "Koefisien Beta3") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x4), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x4", fill = "Koefisien Beta4") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x6), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x5", fill = "Koefisien Beta5") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x7), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x6", fill = "Koefisien Beta6") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x8), color = NA) + # Gunakan 'fill' untuk mengisi warna berdasarkan x1
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x7", fill = "Koefisien Beta7") +
theme_minimal()
# Membuat peta dengan geom_sf
ggplot(data = jabar_merged1) +
geom_sf(aes(fill = x10), color = NA) +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabupaten.kota), size = 3, color = "white") +
labs(title = "Koefisien GWR untuk x8", fill = "Koefisien Beta8") +
theme_minimal()
## Penentuan Kelompok Variabel Signifikan Tiap Kabupaten/Kota
# Ambil hasil koefisien dan standar error
hasil_gwr <- as.data.frame(model_gwr$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 35 variables:
## $ sum.w : num 17.4 17.2 17.3 16.6 16 ...
## $ X.Intercept. : num -1214 -1235 -1252 -1221 -1128 ...
## $ x1 : num 0.00258 0.00252 0.00254 0.00258 0.00261 ...
## $ x2 : num -32.4 -32.8 -33.5 -37.1 -39.5 ...
## $ x3 : num 842 842 834 917 1054 ...
## $ x4 : num 0.0135 0.014 0.0141 0.0177 0.0218 ...
## $ x6 : num 1.6 1.59 1.61 1.44 1.13 ...
## $ x7 : num 8.99 9.26 9.48 9.13 7.85 ...
## $ x8 : num 0.000311 0.000319 0.000315 0.000335 0.000362 ...
## $ x10 : num 0.017 0.016 0.0157 0.0149 0.016 ...
## $ X.Intercept._se : num 518 518 518 514 517 ...
## $ x1_se : num 0.00177 0.00177 0.00177 0.00177 0.00177 ...
## $ x2_se : num 17.6 17.6 17.8 17.6 17.4 ...
## $ x3_se : num 280 282 282 280 280 ...
## $ x4_se : num 0.0118 0.0118 0.0118 0.0118 0.012 ...
## $ x6_se : num 1.09 1.09 1.09 1.09 1.11 ...
## $ x7_se : num 5.25 5.24 5.26 5.2 5.22 ...
## $ x8_se : num 0.000223 0.000223 0.000224 0.000225 0.000228 ...
## $ x10_se : num 0.0115 0.0115 0.0115 0.0115 0.0115 ...
## $ gwr.e : num 175.6 225 10.5 -157.9 -136.4 ...
## $ pred : num 592 382 103 642 281 ...
## $ pred.se : num 86.6 57 53.4 75.2 81 ...
## $ localR2 : num 0.847 0.839 0.838 0.836 0.849 ...
## $ X.Intercept._se_EDF: num 542 542 543 539 542 ...
## $ x1_se_EDF : num 0.00185 0.00185 0.00186 0.00185 0.00186 ...
## $ x2_se_EDF : num 18.5 18.4 18.6 18.4 18.3 ...
## $ x3_se_EDF : num 293 295 296 294 294 ...
## $ x4_se_EDF : num 0.0123 0.0124 0.0124 0.0124 0.0126 ...
## $ x6_se_EDF : num 1.14 1.14 1.14 1.14 1.16 ...
## $ x7_se_EDF : num 5.5 5.49 5.51 5.45 5.47 ...
## $ x8_se_EDF : num 0.000234 0.000234 0.000235 0.000235 0.000239 ...
## $ x10_se_EDF : num 0.0121 0.0121 0.0121 0.012 0.012 ...
## $ pred.se.1 : num 90.7 59.8 55.9 78.7 84.9 ...
## $ coord.x : num -6479679 -6915727 -6822558 -7012851 -7202988 ...
## $ coord.y : num 1.07e+08 1.07e+08 1.07e+08 1.08e+08 1.08e+08 ...
library(dplyr)
# 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
# Hitung p-value manual (z = koefisien / se)
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))),
p_value_x4 = 2 * (1 - pnorm(abs(coef_x4 / x4_se))),
p_value_x6 = 2 * (1 - pnorm(abs(coef_x6 / x6_se))),
p_value_x7 = 2 * (1 - pnorm(abs(coef_x7 / x7_se))),
p_value_x8 = 2 * (1 - pnorm(abs(coef_x8 / x8_se))),
p_value_x10 = 2 * (1 - pnorm(abs(coef_x10 / x10_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),
signifikan_x4 = ifelse(p_value_x4 < 0.05, "x4",NA),
signifikan_x6 = ifelse(p_value_x6 < 0.05, "x6",NA),
signifikan_x7 = ifelse(p_value_x7 < 0.05, "x7",NA),
signifikan_x8 = ifelse(p_value_x8 < 0.05, "x8",NA),
signifikan_x10 = ifelse(p_value_x10 < 0.05, "x10",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, signifikan_x4,
signifikan_x6, signifikan_x7, signifikan_x8, signifikan_x10, 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 x3
## 2 sukabumi x3
## 3 cianjur x3
## 4 bandung x2, x3
## 5 garut x2, x3
## 6 tasikmalaya x2, x3
## 7 ciamis x2, x3
## 8 kuningan x2, x3
## 9 cirebon x2, x3
## 10 majalengka x2, x3
## 11 sumedang x2, x3
## 12 indramayu x2, x3
## 13 subang x2, x3
## 14 purwakarta x3
## 15 karawang x2, x3
## 16 bekasi x3
## 17 bandung barat x2, x3
## 18 pangandaran x2, x3
## 19 kota bogor x3
## 20 kota sukabumi x3
## 21 kota bandung x2, x3
## 22 kota cirebon x2, x3
## 23 kota bekasi x3
## 24 kota depok x3
## 25 kota cimahi x3
## 26 kota tasikmalaya x2, x3
## 27 kota banjar x2, x3
view(tabel_hasil)