UTS FIX SPASIAL
Suhaila
2024-10-17
R Markdown
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:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00R Markdown
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:
Including Plots
You can also embed plots, for example:
## 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 geometriesggplot() +
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 4 tetangga terdekat)
k_neighbors1 <- 4
knn_weights1 <- knearneigh(coords, k = k_neighbors1)
knn_nb1 <- knn2nb(knn_weights1)
knn_listw1 <- nb2listw(knn_nb1, style = "W") # Menggunakan bobot standar
k_neighbors2 <- 5
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.8465, p-value = 0.8014
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.13666482 -0.03846154 0.01345865knn_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 = -0.86937, p-value = 0.8077
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.123347931 -0.038461538 0.009533903# 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 4","KNN 5", "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 4 -0.1366648 -0.03846154 0.013458653 0.8013621
## 2 KNN 5 -0.1233479 -0.03846154 0.009533903 0.8076765
## 3 Queen -0.3344744 -0.03846154 0.019338634 0.9833570
## 4 Rook -0.3344744 -0.03846154 0.019338634 0.9833570#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##Heterogenitas##
bptest(model_ols)##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 16.221, df = 8, p-value = 0.03932##Multikolinearitity##
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######GRW#####
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
)#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# 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##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
nb4 <- knn2nb(knearneigh(coords, k=4))
listw4 <- nb2listw(nb4, style="W")
nb5 <- knn2nb(knearneigh(coords, k=5))
listw5 <- nb2listw(nb5, 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_moran4 <- moran.test(residual_gwr, listw4)
print(knn_moran4)##
## Moran I test under randomisation
##
## data: residual_gwr
## weights: listw4
##
## Moran I statistic standard deviate = 1.1164, p-value = 0.1321
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.09101545 -0.03846154 0.01345011knn_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.009527975knn_moran6 <- moran.test(residual_gwr, queen_listw)
print(knn_moran6)##
## 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.01932658knn_moran7 <- moran.test(residual_gwr, rook_listw)
print(knn_moran7)##
## 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##uji normality##
shapiro.test(residual_gwr)##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.94385, p-value = 0.1516##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##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))
##Plot koefisien regresi lokal untuk variabel x3##
spplot(model_gwr$SDF, "x3", main = "Koefisien GWR untuk x3", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x4##
spplot(model_gwr$SDF, "x4", main = "Koefisien GWR untuk x4", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x6##
spplot(model_gwr$SDF, "x6", main = "Koefisien GWR untuk x6", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x7##
spplot(model_gwr$SDF, "x7", main = "Koefisien GWR untuk x7", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x8##
spplot(model_gwr$SDF, "x8", main = "Koefisien GWR untuk x8", col.regions = terrain.colors(100))
##Plot koefisien regresi lokal untuk variabel x9##
spplot(model_gwr$SDF, "x10", main = "Koefisien GWR untuk x10", col.regions = terrain.colors(100))
##Membuat plot -gglot2##
gwr_df <- as.data.frame(model_gwr$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()
# Plot koefisien lokal untuk x4
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x4), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x4", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x6
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x6), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x6", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x7
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x7), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x7", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x8
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x8), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x8", color = "Koefisien") +
theme_minimal()
# Plot koefisien lokal untuk x10
ggplot(gwr_df, aes(x = data$longitude, y = data$latitude)) +
geom_point(aes(color = x10), size = 2) +
scale_color_viridis_c() +
labs(title = "Koefisien GWR untuk x10", color = "Koefisien") +
theme_minimal()
# 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, x3view(tabel_hasil)