Nama: Aisya Putri Syanurli
NPM: 140610220003
email: aisya22002@mail.unpad.ac.id
Kelas: Spatial Statistics
Dosen pengampu: Dr. I Gede Nyoman Mindra Jaya, M.Si.
Tujuan: Rpubs ini berisikan syntax untuk Pemodelan Faktor-faktor Indeks Kesehatan di Jawa Barat Menggunakan Geographically Weighted Regression (GWR).
Langkah-langkah Analisis GWR:
Loading Library
library(ggsn)
## Loading required package: ggplot2
## Loading required package: grid
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
## Checking rgeos availability: FALSE
library(devtools)
## Loading required package: usethis
library(ggplot2)
library(viridisLite)
library(gridExtra)
library(knitr)
library(gt)
## Warning: package 'gt' was built under R version 4.4.1
library(terra)
## terra 1.7.71
##
## Attaching package: 'terra'
## The following object is masked from 'package:knitr':
##
## spin
## The following object is masked from 'package:ggsn':
##
## north
## The following object is masked from 'package:grid':
##
## depth
library(geodata)
## Warning: package 'geodata' was built under R version 4.4.1
library(leaflet)
library(maps)
library(raster)
## Warning: package 'raster' was built under R version 4.4.1
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:ggsn':
##
## scalebar
library(stars)
## Loading required package: abind
## Loading required package: sf
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(maptools)
##
## Attaching package: 'maptools'
## The following object is masked from 'package:sp':
##
## sp2Mondrian
library(gtools)
library(sp)
library(spdep)
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.4.1
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(ggsn)
library(splancs)
## Warning: package 'splancs' was built under R version 4.4.1
##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysis
##
## Attaching package: 'splancs'
## The following object is masked from 'package:raster':
##
## zoom
## The following objects are masked from 'package:terra':
##
## as.points, is.points, zoom
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ tidyr::extract() masks raster::extract(), terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::map() masks maps::map()
## ✖ dplyr::select() masks raster::select()
## ✖ dplyr::tribble() masks tidyr::tribble(), tibble::tribble(), splancs::tribble()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
require(mgcv)
## Loading required package: mgcv
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## The following object is masked from 'package:raster':
##
## getData
##
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'mgcv'
##
## The following object is masked from 'package:gtools':
##
## scat
require(dplyr)
require(DAAG)
## Loading required package: DAAG
## Warning: package 'DAAG' was built under R version 4.4.1
##
## Attaching package: 'DAAG'
##
## The following object is masked from 'package:maps':
##
## ozone
require(reshape2)
## Loading required package: reshape2
##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths
library(gstat)
library(oce)
## Warning: package 'oce' was built under R version 4.4.1
## Loading required package: gsw
## Warning: package 'gsw' was built under R version 4.4.1
##
## Attaching package: 'oce'
##
## The following object is masked from 'package:terra':
##
## rescale
library(dplyr)
library(sf)
library(spatialreg)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## 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
membaca data
setwd('C:\\Users\\Aisya\\Downloads\\')
data <-read.csv("Variabel UTS Spasial.csv", header = TRUE, sep = ";")
head(data)
## Kabupaten.Kota Indeks.Kesehatan Umur.Harapan.Hidup
## 1 Bandung 83.090 73.975
## 2 Bandung Barat 81.220 72.775
## 3 Banjar 79.220 71.460
## 4 Bekasi 83.140 74.020
## 5 Bogor 79.460 71.630
## 6 Ciamis 80.145 72.075
## Bayi.Berat.Badan.Lahir.Rendah...BBLR.
## 1 1639
## 2 684
## 3 124
## 4 518
## 5 1927
## 6 553
## Persentase.Penduduk.Tidak.Berobat.Jalan..Mengobati.sendiri.
## 1 78.80
## 2 80.97
## 3 94.01
## 4 77.63
## 5 67.42
## 6 84.39
## Jumlah.Balita.dengan.Gizi.Kurang Jumlah.Desa.Kelurahan.Yang.Memiliki.Apotek
## 1 8605.0 142
## 2 4596.0 64
## 3 721.0 10
## 4 3035.0 92
## 5 12537.0 162
## 6 1323.5 50
## Angka.Penemuan.TBC Angka.Kesakitan.DBD.per.100.000.Penduduk
## 1 5839.0 110.0
## 2 1771.0 74.0
## 3 274.0 60.0
## 4 4867.0 25.0
## 5 12153.0 31.0
## 6 1042.5 90.5
## Persentase.Penduduk.Merokok Persentase.Penduduk.Memiliki.BPJS.Kesehatan..PBI.
## 1 15.260 31.130
## 2 17.320 29.180
## 3 11.170 29.670
## 4 10.130 18.670
## 5 14.100 28.400
## 6 13.325 23.985
## Rumah.Tangga.Sanitasi.Layak Jumlah.HIV.AIDS Jumlah.Penyakit.Kusta
## 1 70.85 40 0.050
## 2 50.96 7 0.350
## 3 88.26 10 0.150
## 4 89.99 226 7.570
## 5 73.47 251 3.780
## 6 79.25 45 1.635
attach(data)
Plot Jawa Barat
Indo_Kec<-readRDS('gadm36_IDN_2_sp.rds')
Jabar<-Indo_Kec[Indo_Kec$NAME_1 == "Jawa Barat",]
Jabar$NAME_2
## [1] "Bandung" "Bandung Barat" "Banjar" "Bekasi"
## [5] "Bogor" "Ciamis" "Cianjur" "Cimahi"
## [9] "Cirebon" "Depok" "Garut" "Indramayu"
## [13] "Karawang" "Kota Bandung" "Kota Bekasi" "Kota Bogor"
## [17] "Kota Cirebon" "Kota Sukabumi" "Kota Tasikmalaya" "Kuningan"
## [21] "Majalengka" "Purwakarta" "Subang" "Sukabumi"
## [25] "Sumedang" "Tasikmalaya" "Waduk Cirata"
Jabar <- Jabar[Jabar$NAME_2 != "Waduk Cirata",]
plot(Jabar)
Dataset ini menggunakan data Indeks kesehatan sebagai variabel dependen dan data jumlah balita dengan gizi kurang, jumlah desa kelurahan yang memiliki apotek, persentase penduduk merokok, dan rumah tangga sanitasi layak sebagai variabel independen.
Jabar$Y0 <- as.matrix(data$Indeks.Kesehatan)
Jabar$X1 <- as.matrix(data$Jumlah.Balita.dengan.Gizi.Kurang)
Jabar$X2 <- as.matrix(data$Jumlah.Desa.Kelurahan.Yang.Memiliki.Apotek)
Jabar$X3 <- as.matrix(data$Persentase.Penduduk.Merokok)
Jabar$X4 <- as.matrix(data$Rumah.Tangga.Sanitasi.Layak)
Jabar$Kab.Kota <- as.matrix(data$Kabupaten.Kota)
dilakukan ekplorasi dengan melihat nilai min, median, mean, dan max dari data
Y0 <- Jabar$Y0
x1 <- Jabar$X1
x2 <- Jabar$X2
x3 <- Jabar$X3
x4 <- Jabar$X4
data1=data.frame(Y0,x1,x2,x3,x4)
summary(data1)
## Y0 x1 x2 x3
## Min. :76.85 Min. : 721 Min. : 10.00 Min. :10.13
## 1st Qu.:79.54 1st Qu.: 1725 1st Qu.: 50.00 1st Qu.:12.59
## Median :81.06 Median : 2956 Median : 79.00 Median :13.23
## Mean :81.14 Mean : 3759 Mean : 78.27 Mean :13.70
## 3rd Qu.:83.11 3rd Qu.: 4475 3rd Qu.: 94.25 3rd Qu.:15.06
## Max. :85.35 Max. :12537 Max. :162.00 Max. :18.76
## x4
## Min. :45.80
## 1st Qu.:62.15
## Median :78.22
## Mean :73.36
## 3rd Qu.:86.25
## Max. :96.21
Sebelum melakukan analisis GWR, maka perlu diuji analisis yang lebih sederhana terlebih dahulu. Lakukan analisis regresi meggunakan metode OLS untuk melihat pengaruh faktor indeks kesehatan untuk wilayah jawa barat secara global. Regresi Global Metode OLS perlu menguji asumsi klasik sebelum menggunakan analisisnya.
Multikolinearitas adalah kondisi di mana terdapat hubungan linear yang kuat antara dua atau lebih variabel independen dalam model regresi. Multikolinearitas dapat dilihat menggunakan nilai VIF (Variance Inflation Factor).
$$ \text{VIF}_i = \frac{1}{1 - R_i^2} $$
Nilai VIF yang > 5 atau > 10 mengindikasikan ada nya multikolinearitas antar variabel
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:DAAG':
##
## vif
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:gtools':
##
## logit
## The following object is masked from 'package:maptools':
##
## pointLabel
jabarlm <- lm(Y0 ~ (x1+x2+x3+x4))
vif(jabarlm)
## x1 x2 x3 x4
## 2.904493 2.591739 1.856470 1.608910
Dalam model regresi yang ideal, diasumsikan bahwa error memiliki varians yang konstan (homoskedastisitas), yang berarti penyebaran error sekitar garis regresi seharusnya sama di semua titik. Jika varians error berubah-ubah atau berbeda-beda pada berbagai tingkat variabel independen, maka terjadi heteroskedastisitas.
Hipotesis:
H₀: Tidak ada heteroskedastisitas dalam model, atau varians error adalah konstan di seluruh rentang variabel independen (homoskedastisitas).
H₁: Terdapat heteroskedastisitas dalam model, atau varians error tidak konstan di seluruh rentang variabel independen.
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
##
## time<-
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
bptest(jabarlm)
##
## studentized Breusch-Pagan test
##
## data: jabarlm
## BP = 8.2735, df = 4, p-value = 0.08206
Normalitas residual adalah asumsi bahwa residual (galat atau error) dalam model regresi linier berdistribusi normal. Jika residual tidak berdistribusi normal, maka estimasi koefisien tetap tidak bias, namun hasil uji statistik menjadi tidak valid, yang dapat menyebabkan kesimpulan yang salah.
Hipotesis:
H₀: Residual dalam model regresi berdistribusi normal.
H₁: Residual dalam model regresi tidak berdistribusi normal.
resids <- residuals(jabarlm)
shapiro.test(resids)
##
## Shapiro-Wilk normality test
##
## data: resids
## W = 0.97781, p-value = 0.8243
Autokorelasi adalah kondisi di mana residual dari suatu observasi memiliki korelasi dengan residual dari observasi lainnya. Autokorelasi biasanya terjadi dalam data runtun waktu (time series) di mana nilai residual memiliki pola berurutan. Asumsi non-autokorelasi (independensi residual) menyatakan bahwa residual dari observasi satu tidak berkorelasi dengan residual dari observasi lainnya.
Jika autokorelasi hadir, maka model regresi mungkin tidak tepat dan estimasi koefisien menjadi bias, yang mengganggu validitas inferensi statistik.
Hipotesis:
H₀: Tidak ada autokorelasi dalam residual (residual independen)
H₁: Terdapat autokorelasi dalam residual.
durbinWatsonTest(jabarlm)
## lag Autocorrelation D-W Statistic p-value
## 1 0.04869874 1.650291 0.368
## Alternative hypothesis: rho != 0
jabarlm <- lm(Y0 ~ (x1+x2+x3+x4))
summary(jabarlm)
##
## Call:
## lm(formula = Y0 ~ (x1 + x2 + x3 + x4))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6270 -1.3180 0.0582 1.6168 3.3910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.456e+01 6.359e+00 13.296 1.07e-11 ***
## x1 5.465e-05 2.656e-04 0.206 0.839
## x2 -5.238e-03 1.713e-02 -0.306 0.763
## x3 -2.773e-01 3.142e-01 -0.883 0.387
## x4 8.083e-03 3.807e-02 0.212 0.834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.278 on 21 degrees of freedom
## Multiple R-squared: 0.08228, Adjusted R-squared: -0.09253
## F-statistic: 0.4707 on 4 and 21 DF, p-value: 0.7566
Berdasarkan hasil analisis menggunakan analisis regresi global menggunakan metode OLS didapatkan koefisien variabel x1, x2, x3, dan x4 tidak signifikan terhadap model. Nilai koefisien determinasi yang dihasilkan adalah 8.2%, artinya variabel faktor-faktor dari indeks kesehatan belum bagus dalam menjelaskan variasi dari model indeks kesehatan tersebut.n
par(mfrow=c(2,2))
plot(jabarlm)
Autokorelasi spasial adalah kecenderungan variabel di lokasi berdekatan untuk memiliki nilai yang lebih mirip dibandingkan dengan variabel di lokasi yang lebih jauh. Uji autokorelasi spasial dapat menggunakan metode Moran’s I. Moran’s I dilakukan dengan melihat kedekatan lokasi suatu pengamatan dengan lokasi pengamatan lain memiliki pengaruh satu sama lain.
Hipotesis:
H₀: (tidak ada autokorelasi antar ruang/lokasi)
H₁: (terdapat autokorelasi antar ruang/lokasi
Jabar$id<-c(1:26)
data$id<-c(1:26)
data$x<-coordinates(Jabar)[,1]
data$Y<-coordinates(Jabar)[,2]
Jabar_sf<-st_as_sf(Jabar)
Jabar_merged <- Jabar_sf %>%
left_join(data, by = "id")
#Plot
ggplot() +
geom_sf(data=Jabar_merged, aes(fill = Indeks.Kesehatan),color=NA) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
theme(legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank())+ # Remove Y0-axis labels
labs(title = "",
fill = "data")
membuat peta interval
summary(data$Indeks.Kesehatan)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 76.85 79.54 81.06 81.14 83.11 85.35
breaks <- c(-Inf, 79.61, 80.96, 83.11,Inf)
# Define labels for each interval
labels <- c("Very Low", "Low", "High",
"Very High")
# Create a new column with discretized data
Jabar_merged$data_Discrete <-
cut(Jabar_merged$Indeks.Kesehatan, breaks =
breaks, labels = labels, right = TRUE)
ggplot() +
geom_sf(data=Jabar_merged, aes(fill =
data_Discrete),color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow",
"Low" = "orange",
"High" = "red",
"Very High" = "red3"))+
labs(fill = "data")+theme(legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank())+ # Remove Y0-axis labels
labs(title = "",
fill = "data")
ID <- c(1:26)
CoordK<-coordinates(Jabar)
row.names(Jabar) = as.character(1:26)
W <- poly2nb(Jabar, row.names=ID,
queen=FALSE) ; W
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 102
## Percentage nonzero weights: 15.08876
## Average number of links: 3.923077
WB <- nb2mat(W, style='B', zero.policy = TRUE); WB #menyajikan dalam bentuk matrix biner "B"
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## 1 0 1 0 0 0 0 1 1 0 0 1 0 0 1
## 2 1 0 0 0 0 0 1 1 0 0 0 0 0 1
## 3 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 4 0 0 0 0 1 0 0 0 0 0 0 0 1 0
## 5 0 0 0 1 0 0 1 0 0 1 0 0 1 0
## 6 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## 7 1 1 0 0 1 0 0 0 0 0 1 0 0 0
## 8 1 1 0 0 0 0 0 0 0 0 0 0 0 1
## 9 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## 10 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 11 1 0 0 0 0 0 1 0 0 0 0 0 0 0
## 12 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 13 0 0 0 1 1 0 0 0 0 0 0 0 0 0
## 14 1 1 0 0 0 0 0 1 0 0 0 0 0 0
## 15 0 0 0 1 1 0 0 0 0 1 0 0 0 0
## 16 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 17 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 18 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 19 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 20 0 0 0 0 0 1 0 0 1 0 0 0 0 0
## 21 0 0 0 0 0 1 0 0 1 0 0 1 0 0
## 22 0 1 0 0 1 0 1 0 0 0 0 0 1 0
## 23 1 1 0 0 0 0 0 0 0 0 0 1 1 0
## 24 0 0 0 0 1 0 1 0 0 0 0 0 0 0
## 25 1 0 0 0 0 0 0 0 0 0 1 1 0 0
## 26 0 0 0 0 0 1 0 0 0 0 1 0 0 0
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## 1 0 0 0 0 0 0 0 0 1 0 1 0
## 2 0 0 0 0 0 0 0 1 1 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0 0 0 0 0
## 5 1 1 0 0 0 0 0 1 0 1 0 0
## 6 0 0 0 0 1 1 1 0 0 0 0 1
## 7 0 0 0 0 0 0 0 1 0 1 0 0
## 8 0 0 0 0 0 0 0 0 0 0 0 0
## 9 0 0 1 0 0 1 1 0 0 0 0 0
## 10 1 0 0 0 0 0 0 0 0 0 0 0
## 11 0 0 0 0 0 0 0 0 0 0 1 1
## 12 0 0 0 0 0 0 1 0 1 0 1 0
## 13 0 0 0 0 0 0 0 1 1 0 0 0
## 14 0 0 0 0 0 0 0 0 0 0 0 0
## 15 0 0 0 0 0 0 0 0 0 0 0 0
## 16 0 0 0 0 0 0 0 0 0 0 0 0
## 17 0 0 0 0 0 0 0 0 0 0 0 0
## 18 0 0 0 0 0 0 0 0 0 1 0 0
## 19 0 0 0 0 0 0 0 0 0 0 0 1
## 20 0 0 0 0 0 0 1 0 0 0 0 0
## 21 0 0 0 0 0 1 0 0 0 0 1 1
## 22 0 0 0 0 0 0 0 0 1 0 0 0
## 23 0 0 0 0 0 0 0 1 0 0 1 0
## 24 0 0 0 1 0 0 0 0 0 0 0 0
## 25 0 0 0 0 0 0 1 0 1 0 0 1
## 26 0 0 0 0 1 0 1 0 0 0 1 0
## attr(,"call")
## nb2mat(neighbours = W, style = "B", zero.policy = TRUE)
WL_Rook <- nb2listw(W); WL_Rook #List neighbours
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 102
## Percentage nonzero weights: 15.08876
## Average number of links: 3.923077
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 26 676 26 16.19885 115.1439
plot(Jabar, axes=T, col="white")
text(CoordK[,1], CoordK[,2],
row.names(Jabar), col="black", cex=0.8,
pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19,
cex=0.5,col="blue")
#global moran Rook
MI_Rook <- moran.test(data$Indeks.Kesehatan, WL_Rook)
moran.plot(data$Indeks.Kesehatan, WL_Rook)
#local moran Rook
locm<-localmoran(data$Indeks.Kesehatan, WL_Rook)
head(locm)
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.10375296 -0.0331742566 0.0893482348 0.4580858 0.64689083
## 2 0.00759965 -0.0000510429 0.0001750966 0.5781785 0.56314363
## 3 0.42057451 -0.0324052135 0.8152330072 0.5016928 0.61588361
## 4 0.33930659 -0.0349005870 0.2675890365 0.7233991 0.46943473
## 5 -0.20046101 -0.0248237078 0.0557276624 -0.7440144 0.45686776
## 6 0.32656225 -0.0087335537 0.0375148745 1.7311168 0.08343093
Perbedaan garisnya lebih banyak
W <- poly2nb(Jabar, row.names=ID,
queen=TRUE) #Mendapatkan W
WB <- nb2mat(W, style='B', zero.policy =
TRUE) #menyajikan dalam bentuk matrix biner "B"
WB[1:5,1:5]
## [,1] [,2] [,3] [,4] [,5]
## 1 0 1 0 0 0
## 2 1 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 1
## 5 0 0 0 1 0
WL_Queen<-nb2listw(W); WL_Queen
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 102
## Percentage nonzero weights: 15.08876
## Average number of links: 3.923077
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 26 676 26 16.19885 115.1439
plot(Jabar, axes=T, col="white")
text(CoordK[,1], CoordK[,2],
row.names(Jabar), col="black", cex=0.8,
pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19,
cex=0.5,col="blue")
#global moran Queen
MI_Queen = moran.test(data$Indeks.Kesehatan, WL_Queen)
moran.plot(data$Indeks.Kesehatan, WL_Queen)
#local moran Queen
locm<-localmoran(data$Indeks.Kesehatan, WL_Queen)
head(locm)
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.10375296 -0.0331742566 0.0893482348 0.4580858 0.64689083
## 2 0.00759965 -0.0000510429 0.0001750966 0.5781785 0.56314363
## 3 0.42057451 -0.0324052135 0.8152330072 0.5016928 0.61588361
## 4 0.33930659 -0.0349005870 0.2675890365 0.7233991 0.46943473
## 5 -0.20046101 -0.0248237078 0.0557276624 -0.7440144 0.45686776
## 6 0.32656225 -0.0087335537 0.0375148745 1.7311168 0.08343093
#3. Berdasarkan Jarak
#K-nearesr neigbhors k = 3
as.matrix(CoordK)
## [,1] [,2]
## 487 107.6108 -7.099969
## 476 107.4150 -6.897056
## 496 108.5665 -7.376302
## 497 107.1207 -6.215149
## 498 106.7687 -6.561184
## 499 108.4661 -7.434347
## 500 107.1578 -7.133713
## 501 107.5436 -6.886495
## 502 108.5513 -6.745512
## 477 106.8168 -6.396101
## 478 107.7889 -7.359586
## 479 108.1687 -6.448640
## 480 107.3539 -6.252003
## 481 107.6366 -6.919135
## 482 106.9757 -6.280230
## 483 106.7996 -6.593453
## 484 108.5535 -6.741886
## 485 106.9243 -6.939871
## 486 108.2194 -7.360251
## 488 108.5603 -7.004233
## 489 108.2578 -6.815865
## 490 107.4322 -6.595016
## 491 107.7322 -6.484194
## 492 106.7101 -7.074623
## 493 107.9808 -6.825066
## 494 108.1413 -7.496892
W <- knn2nb(knearneigh(CoordK, k =
3), row.names = ID); W
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 78
## Percentage nonzero weights: 11.53846
## Average number of links: 3
## Non-symmetric neighbours list
WB <- nb2mat(W, style='B', zero.policy
= TRUE) #menyajikan dalam bentuk matrix biner "B"
WB[1:5,1:5]
## [,1] [,2] [,3] [,4] [,5]
## 1 0 1 0 0 0
## 2 1 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
WL_k3 <- nb2listw(W); WL_k3 #untuk Morans's I
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 78
## Percentage nonzero weights: 11.53846
## Average number of links: 3
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 26 676 26 15.55556 107.7778
plot(Jabar, axes=T, col="white")
text(CoordK[,1], CoordK[,2],
row.names(Jabar), col="black", cex=0.8,
pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19,
cex=0.5,col="blue")
#global moran K-nearesr neigbhors k = 3
MI_k3 <- moran.test(data$Indeks.Kesehatan, WL_k3); MI_k3
##
## Moran I test under randomisation
##
## data: data$Indeks.Kesehatan
## weights: WL_k3
##
## Moran I statistic standard deviate = 2.2216, p-value = 0.01315
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27559679 -0.04000000 0.02018025
moran.plot(data$Indeks.Kesehatan, WL_k3)
#local moran k = 3
locm<-localmoran(data$Indeks.Kesehatan, WL_k3)
head(locm)
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.83361751 -0.0331742566 0.2548079288 1.717150 0.085951715
## 2 0.04311754 -0.0000510429 0.0004054868 2.143776 0.032050836
## 3 -0.11287130 -0.0324052135 0.2490989744 -0.161223 0.871917787
## 4 1.07210349 -0.0349005870 0.2675890365 2.140006 0.032354253
## 5 -1.19045581 -0.0248237078 0.1923150701 -2.657999 0.007860609
## 6 0.46575870 -0.0087335537 0.0687772700 1.809284 0.070406948
#K-nearesr neigbhors k = 5
W <- knn2nb(knearneigh(CoordK, k =
5), row.names = ID); W
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 130
## Percentage nonzero weights: 19.23077
## Average number of links: 5
## Non-symmetric neighbours list
WB <- nb2mat(W, style='B', zero.policy
= TRUE) #menyajikan dalam bentuk matrix biner "B"
WB[1:5,1:5]
## [,1] [,2] [,3] [,4] [,5]
## 1 0 1 0 0 0
## 2 1 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 1
## 5 0 0 0 1 0
WL_k5 <- nb2listw(W); WL_k5 #untuk Morans's I
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 26
## Number of nonzero links: 130
## Percentage nonzero weights: 19.23077
## Average number of links: 5
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 26 676 26 9.76 105.6
plot(Jabar, axes=T, col="white")
text(CoordK[,1], CoordK[,2],
row.names(Jabar), col="black", cex=0.8,
pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19,
cex=0.5,col="blue")
#global moran K-nearesr neigbhors k = 5
MI_k5 <- moran.test(data$Indeks.Kesehatan, WL_k5); MI_k5
##
## Moran I test under randomisation
##
## data: data$Indeks.Kesehatan
## weights: WL_k5
##
## Moran I statistic standard deviate = 1.1473, p-value = 0.1256
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.08321328 -0.04000000 0.01153410
moran.plot(data$Indeks.Kesehatan, WL_k5)
#local moran k = 5
locm<-localmoran(data$Indeks.Kesehatan, WL_k5)
head(locm)
## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.099857102 -0.0331742566 0.1389861430 0.3568357 0.721214823
## 2 0.008195793 -0.0000510429 0.0002211746 0.5545232 0.579220844
## 3 0.329607963 -0.0324052135 0.1358721679 0.9821077 0.326046782
## 4 0.324446696 -0.0349005870 0.1459576563 0.9405912 0.346914378
## 5 -0.873706618 -0.0248237078 0.1048991291 -2.6209688 0.008768029
## 6 0.326562255 -0.0087335537 0.0375148745 1.7311168 0.083430932
MI_Rook$estimate
## Moran I statistic Expectation Variance
## 0.28506018 -0.04000000 0.02071532
MI_Queen$estimate
## Moran I statistic Expectation Variance
## 0.28506018 -0.04000000 0.02071532
MI_k3$estimate
## Moran I statistic Expectation Variance
## 0.27559679 -0.04000000 0.02018025
MI_k5$estimate
## Moran I statistic Expectation Variance
## 0.08321328 -0.04000000 0.01153410
cari yang moran I statistik yang paling besar
moran_index<-data.frame(
"Weight matrix"=c("Rook Contiguity", "Queen Contiguity", "KNN (k=3)", "KNN (k=5)"),
"Moran's I"=c(MI_Rook$estimate["Moran I statistic"], MI_Queen$estimate["Moran I statistic"], MI_k3$estimate["Moran I statistic"], MI_k5$estimate["Moran I statistic"]),
"P-value"=c(MI_Rook$p.value, MI_Queen$p.value, MI_k3$p.value, MI_k5$p.value)); moran_index
## Weight.matrix Moran.s.I P.value
## 1 Rook Contiguity 0.28506018 0.01195760
## 2 Queen Contiguity 0.28506018 0.01195760
## 3 KNN (k=3) 0.27559679 0.01315457
## 4 KNN (k=5) 0.08321328 0.12563505
W_optimum <- WL_Queen
Pilih W dari Moran’s I yang signifikan dan terbesar. Dalam kasus ini adalah rook contiguity dan queen contiguity. Sehingga dapat disimpulkan bahwa data memiliki autokorelasi spasial.
library(spgwr)
## Warning: package 'spgwr' was built under R version 4.4.1
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
adapt_gaussian <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar,
coords = coo, adapt = T)
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = coo, adapt =
## T): data is Spatial* object, ignoring coords argument
## Adaptive q: 0.381966 CV score: 170.729
## Adaptive q: 0.618034 CV score: 172.887
## Adaptive q: 0.236068 CV score: 172.8794
## Adaptive q: 0.4268928 CV score: 169.4459
## Adaptive q: 0.4999023 CV score: 170.2355
## Adaptive q: 0.4472012 CV score: 169.3491
## Adaptive q: 0.4451121 CV score: 169.3569
## Adaptive q: 0.4673312 CV score: 169.4376
## Adaptive q: 0.4512578 CV score: 169.3355
## Adaptive q: 0.4573973 CV score: 169.3189
## Adaptive q: 0.4611917 CV score: 169.3111
## Adaptive q: 0.4635368 CV score: 169.3537
## Adaptive q: 0.4597424 CV score: 169.3139
## Adaptive q: 0.4620875 CV score: 169.3223
## Adaptive q: 0.4606381 CV score: 169.3121
## Adaptive q: 0.4615339 CV score: 169.3105
## Adaptive q: 0.4617453 CV score: 169.3149
## Adaptive q: 0.4614032 CV score: 169.3107
## Adaptive q: 0.4616146 CV score: 169.3121
## Adaptive q: 0.4614839 CV score: 169.3106
## Adaptive q: 0.4615339 CV score: 169.3105
adapt_gaussian
## [1] 0.4615339
gwr.model1 = gwr(Y0 ~ X1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
adapt=adapt_gaussian,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model1
## Warning in gwr(Y0 ~ X1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ X1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], adapt = adapt_gaussian, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.4615339 (about 11 of 26 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 7.9506e+01 8.3195e+01 8.7902e+01 9.2346e+01 9.3701e+01
## X1 -8.4330e-05 7.7633e-06 1.0531e-04 1.3347e-04 1.7254e-04
## x2 -2.1113e-02 -1.6694e-02 -6.2342e-03 2.4634e-03 8.2923e-03
## x3 -7.0297e-01 -6.4822e-01 -4.1455e-01 -1.4718e-01 1.5859e-02
## x4 -2.7339e-02 -1.2292e-02 -7.4870e-03 -1.6419e-03 9.8735e-03
## Global
## X.Intercept. 84.5552
## X1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 9.422341
## Effective degrees of freedom (residual: 2traceS - traceS'S): 16.57766
## Sigma (residual: 2traceS - traceS'S): 1.998207
## Effective number of parameters (model: traceS): 7.716878
## Effective degrees of freedom (model: traceS): 18.28312
## Sigma (model: traceS): 1.902729
## Sigma (ML): 1.595568
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 125.918
## AIC (GWR p. 96, eq. 4.22): 105.7976
## Residual sum of squares: 66.19178
## Quasi-global R2: 0.4426452
results1 <-as.data.frame(gwr.model1$SDF)
head(results1)
## sum.w X.Intercept. X1 x2 x3 x4
## 1 15.39218 86.12390 8.323564e-05 -0.002674802 -0.330945477 -0.006355016
## 2 16.17133 89.67487 1.313558e-04 -0.009942713 -0.510627029 -0.012460061
## 3 14.56258 80.17200 -6.502295e-05 0.006632108 -0.007770284 0.005986146
## 4 14.91016 92.38294 1.451677e-04 -0.018727817 -0.651593937 -0.011787012
## 5 14.57546 92.23680 1.322143e-04 -0.019624322 -0.662271117 -0.005983475
## 6 14.91581 80.16839 -6.257390e-05 0.006545607 -0.012301396 0.006984880
## X.Intercept._se X1_se x2_se x3_se x4_se gwr.e
## 1 5.468670 0.0002359812 0.01531514 0.2700406 0.03300336 2.1301621
## 2 5.638181 0.0002276313 0.01523083 0.2747881 0.03380755 1.0567737
## 3 5.731380 0.0002616676 0.01592777 0.2863028 0.03489529 -1.4129831
## 4 6.064310 0.0002323723 0.01567795 0.2936172 0.03569832 -0.2992089
## 5 6.114433 0.0002332142 0.01583684 0.2921314 0.03644730 -1.4776068
## 6 5.714456 0.0002584592 0.01580563 0.2850959 0.03467637 -0.6574872
## pred pred.se localR2 X.Intercept._se_EDF X1_se_EDF x2_se_EDF
## 1 80.95984 0.8232101 0.3783692 5.743086 0.0002478226 0.01608365
## 2 80.16323 0.8613406 0.4427207 5.921103 0.0002390538 0.01599511
## 3 80.63298 0.9859878 0.2796241 6.018978 0.0002747980 0.01672702
## 4 83.43921 0.9184406 0.5347289 6.368615 0.0002440326 0.01646467
## 5 80.93761 1.3968399 0.5391102 6.421253 0.0002449168 0.01663153
## 6 80.80249 0.5774716 0.2862436 6.001205 0.0002714286 0.01659875
## x3_se_EDF x4_se_EDF pred.se.1
## 1 0.2835911 0.03465945 0.8645185
## 2 0.2885768 0.03550400 0.9045624
## 3 0.3006693 0.03664633 1.0354642
## 4 0.3083507 0.03748965 0.9645276
## 5 0.3067905 0.03827621 1.4669328
## 6 0.2994019 0.03641642 0.6064489
Jabar$bwadapt <- gwr.model1$bandwidth
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
tm_shape(Jabar, unit = "mi") +
tm_polygons(col = "bwadapt", style = "quantile",palette = "Oranges",
border.alpha = 0, title = "") +
tm_scale_bar(breaks = c(0, 1, 2), size = 1, position = c("right", "bottom")) +
tm_borders(col = "black", lwd = 0.5)+
tm_layout(main.title = "GWR bandwidth", main.title.size = 0.95,
frame = FALSE, legend.text.size = 1)
## Warning: The argument size of tm_scale_bar is deprecated. It has been renamed
## to text.size
## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).
adapt_bisquare <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar,
coords = CoordK[c(1:27),c(1,2)], gweight = gwr.bisquare, adapt = T)
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords =
## CoordK[c(1:27), : data is Spatial* object, ignoring coords argument
## Adaptive q: 0.381966 CV score: 180.5357
## Adaptive q: 0.618034 CV score: 183.3186
## Adaptive q: 0.236068 CV score: 1903.262
## Adaptive q: 0.4998095 CV score: 160.0533
## Adaptive q: 0.4962454 CV score: 159.7837
## Adaptive q: 0.4806981 CV score: 159.2194
## Adaptive q: 0.4796629 CV score: 159.2273
## Adaptive q: 0.4816221 CV score: 159.218
## Adaptive q: 0.4814047 CV score: 159.2179
## Adaptive q: 0.4814454 CV score: 159.2179
## Adaptive q: 0.481364 CV score: 159.2179
## Adaptive q: 0.4814047 CV score: 159.2179
adapt_bisquare
## [1] 0.4814047
gwr.model2 = gwr(Y0 ~ x1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
adapt=adapt_bisquare,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model2
## Warning in gwr(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], adapt = adapt_bisquare, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.4814047 (about 12 of 26 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 7.9960e+01 8.3223e+01 8.7815e+01 9.2088e+01 9.3538e+01
## x1 -6.6549e-05 8.7804e-06 1.0489e-04 1.3159e-04 1.7088e-04
## x2 -2.0505e-02 -1.6505e-02 -6.1988e-03 2.4001e-03 7.1747e-03
## x3 -6.8734e-01 -6.3960e-01 -4.1087e-01 -1.5024e-01 -9.9036e-03
## x4 -2.6635e-02 -1.2167e-02 -6.8638e-03 -1.5317e-03 9.1260e-03
## Global
## X.Intercept. 84.5552
## x1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 9.242917
## Effective degrees of freedom (residual: 2traceS - traceS'S): 16.75708
## Sigma (residual: 2traceS - traceS'S): 2.013418
## Effective number of parameters (model: traceS): 7.589667
## Effective degrees of freedom (model: traceS): 18.41033
## Sigma (model: traceS): 1.92089
## Sigma (ML): 1.616392
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 125.9734
## AIC (GWR p. 96, eq. 4.22): 106.3447
## Residual sum of squares: 67.93077
## Quasi-global R2: 0.4280024
results2 <-as.data.frame(gwr.model2$SDF)
head(results2)
## sum.w X.Intercept. x1 x2 x3 x4
## 1 16.81316 85.79908 7.590640e-05 -0.003097754 -0.31788645 -0.003839219
## 2 16.18879 89.66138 1.311494e-04 -0.009932748 -0.51001550 -0.012402658
## 3 15.16154 80.54643 -5.110301e-05 0.005782165 -0.03091073 0.005728613
## 4 15.14246 92.16028 1.436003e-04 -0.018346083 -0.64130438 -0.011200889
## 5 14.75787 92.05194 1.312776e-04 -0.019297844 -0.65312918 -0.005625825
## 6 15.18066 80.34112 -5.627498e-05 0.006151013 -0.02272923 0.006824467
## X.Intercept._se x1_se x2_se x3_se x4_se gwr.e
## 1 5.471349 0.0002332370 0.01511982 0.2701785 0.03293079 2.2005776
## 2 5.690384 0.0002297756 0.01537203 0.2773493 0.03412086 1.0570645
## 3 5.724107 0.0002575103 0.01582829 0.2850869 0.03480485 -1.5077375
## 4 6.080134 0.0002340762 0.01575631 0.2946806 0.03583313 -0.2638851
## 5 6.135414 0.0002350099 0.01592934 0.2935335 0.03659125 -1.4890650
## 6 5.741346 0.0002580840 0.01584935 0.2860869 0.03482397 -0.6671645
## pred pred.se localR2 X.Intercept._se_EDF x1_se_EDF x2_se_EDF
## 1 80.88942 0.8083068 0.3742587 5.734903 0.0002444720 0.01584814
## 2 80.16294 0.8694725 0.4298960 5.964488 0.0002408439 0.01611250
## 3 80.72774 0.9814482 0.2733266 5.999836 0.0002699145 0.01659073
## 4 83.40389 0.9215723 0.5206666 6.373013 0.0002453516 0.01651529
## 5 80.94906 1.4057199 0.5255469 6.430956 0.0002463302 0.01669665
## 6 80.81216 0.5790812 0.2724018 6.017905 0.0002705158 0.01661281
## x3_se_EDF x4_se_EDF pred.se.1
## 1 0.2831930 0.03451706 0.8472428
## 2 0.2907092 0.03576445 0.9113547
## 3 0.2988195 0.03648139 1.0287243
## 4 0.3088753 0.03755920 0.9659642
## 5 0.3076730 0.03835384 1.4734331
## 6 0.2998676 0.03650143 0.6069754
adapt_tricube <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar,
coords = CoordK[c(1:27),c(1,2)], gweight = gwr.tricube, adapt = T)
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords =
## CoordK[c(1:27), : data is Spatial* object, ignoring coords argument
## Adaptive q: 0.381966 CV score: 184.6867
## Adaptive q: 0.618034 CV score: 176.9737
## Adaptive q: 0.763932 CV score: 187.0415
## Adaptive q: 0.5613693 CV score: 162.5054
## Adaptive q: 0.5101762 CV score: 159.8907
## Adaptive q: 0.5222876 CV score: 159.8907
## Adaptive q: 0.5162319 CV score: 159.8907
## Adaptive q: 0.5192598 CV score: 159.8907
## Adaptive q: 0.5181032 CV score: 159.8907
## Adaptive q: 0.5173884 CV score: 159.8907
## Adaptive q: 0.5169467 CV score: 159.8907
## Adaptive q: 0.517525 CV score: 159.8907
## Adaptive q: 0.5175657 CV score: 159.8907
## Adaptive q: 0.5174728 CV score: 159.8907
## Adaptive q: 0.517525 CV score: 159.8907
adapt_tricube
## [1] 0.517525
gwr.model3 = gwr(Y0 ~ x1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
adapt=adapt_tricube,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model3
## Warning in gwr(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], adapt = adapt_tricube, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.517525 (about 13 of 26 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 8.0333e+01 8.3248e+01 8.7738e+01 9.1823e+01 9.3389e+01
## x1 -5.3493e-05 9.7140e-06 1.0451e-04 1.3081e-04 1.6934e-04
## x2 -1.9957e-02 -1.6389e-02 -6.1675e-03 2.3414e-03 6.2794e-03
## x3 -6.7308e-01 -6.3172e-01 -4.0759e-01 -1.5217e-01 -3.1127e-02
## x4 -2.5991e-02 -1.2091e-02 -6.3428e-03 -1.4279e-03 8.5847e-03
## Global
## X.Intercept. 84.5552
## x1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 9.093456
## Effective degrees of freedom (residual: 2traceS - traceS'S): 16.90654
## Sigma (residual: 2traceS - traceS'S): 2.026358
## Effective number of parameters (model: traceS): 7.4847
## Effective degrees of freedom (model: traceS): 18.5153
## Sigma (model: traceS): 1.936324
## Sigma (ML): 1.634018
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 126.0339
## AIC (GWR p. 96, eq. 4.22): 106.8037
## Residual sum of squares: 69.42038
## Quasi-global R2: 0.4154594
results3 <-as.data.frame(gwr.model3$SDF)
head(results3)
## sum.w X.Intercept. x1 x2 x3 x4
## 1 17.93778 85.58080 7.142483e-05 -0.003423853 -0.30966784 -0.002003223
## 2 16.20501 89.64885 1.309578e-04 -0.009923493 -0.50944790 -0.012349388
## 3 15.69328 80.85634 -4.004110e-05 0.005058162 -0.05006499 0.005576933
## 4 15.35523 91.95775 1.421251e-04 -0.017999770 -0.63192980 -0.010667344
## 5 14.92605 91.88285 1.303824e-04 -0.018998999 -0.64477295 -0.005297206
## 6 15.42263 80.49395 -5.080426e-05 0.005798007 -0.03195786 0.006695730
## X.Intercept._se x1_se x2_se x3_se x4_se gwr.e
## 1 5.485004 0.0002322381 0.01503875 0.2708824 0.03296089 2.2482317
## 2 5.734602 0.0002315954 0.01549165 0.2795205 0.03438626 1.0573342
## 3 5.721530 0.0002544370 0.01575895 0.2843651 0.03474363 -1.5910462
## 4 6.091558 0.0002354978 0.01581929 0.2955027 0.03593711 -0.2317166
## 5 6.151087 0.0002365082 0.01600390 0.2946469 0.03670141 -1.4991311
## 6 5.763478 0.0002577043 0.01588371 0.2869004 0.03494283 -0.6764101
## pred pred.se localR2 X.Intercept._se_EDF x1_se_EDF x2_se_EDF
## 1 80.84177 0.8014872 0.3698612 5.740039 0.0002430365 0.01573800
## 2 80.16267 0.8763754 0.4191723 6.001243 0.0002423639 0.01621196
## 3 80.81105 0.9777644 0.2687303 5.987563 0.0002662675 0.01649169
## 4 83.37172 0.9239632 0.5082954 6.374796 0.0002464477 0.01655483
## 5 80.95913 1.4129904 0.5136375 6.437093 0.0002475051 0.01674803
## 6 80.82141 0.5803373 0.2613288 6.031462 0.0002696868 0.01662225
## x3_se_EDF x4_se_EDF pred.se.1
## 1 0.2834776 0.03449347 0.8387539
## 2 0.2925173 0.03598512 0.9171241
## 3 0.2975872 0.03635911 1.0232274
## 4 0.3092426 0.03760808 0.9669247
## 5 0.3083471 0.03840791 1.4786901
## 6 0.3002404 0.03656756 0.6073211
library(GWmodel)
## Warning: package 'GWmodel' was built under R version 4.4.1
## Loading required package: robustbase
##
## Attaching package: 'robustbase'
## The following object is masked from 'package:DAAG':
##
## milk
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
adapt_exponential <- bw.gwr(Y0 ~ x1+x2+x3+x4,
data=Jabar,
approach="CV",
kernel="exponential",
adaptive=T)
## Adaptive bandwidth: 23 CV score: 174.1882
## Adaptive bandwidth: 22 CV score: 173.9999
## Adaptive bandwidth: 20 CV score: 173.7872
## Adaptive bandwidth: 20 CV score: 173.7872
adapt_exponential
## [1] 20
gwr.model4 <- gwr.basic(Y0 ~ x1+x2+x3+x4,
data=Jabar,
bw=adapt_exponential,
kernel="exponential") ; gwr.model4
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2024-11-28 23:46:20.339195
## Call:
## gwr.basic(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, bw = adapt_exponential,
## kernel = "exponential")
##
## Dependent (y) variable: Y0
## Independent variables: x1 x2 x3 x4
## Number of data points: 26
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6270 -1.3180 0.0582 1.6168 3.3910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.456e+01 6.359e+00 13.296 1.07e-11 ***
## x1 5.465e-05 2.656e-04 0.206 0.839
## x2 -5.238e-03 1.713e-02 -0.306 0.763
## x3 -2.773e-01 3.142e-01 -0.883 0.387
## x4 8.083e-03 3.807e-02 0.212 0.834
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 2.278 on 21 degrees of freedom
## Multiple R-squared: 0.08228
## Adjusted R-squared: -0.09253
## F-statistic: 0.4707 on 4 and 21 DF, p-value: 0.7566
## ***Extra Diagnostic information
## Residual sum of squares: 108.9892
## Sigma(hat): 2.131013
## AIC: 123.0468
## AICc: 127.4678
## BIC: 124.1439
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: exponential
## Fixed bandwidth: 20
## Regression points: the same locations as observations are used.
## Distance metric: Euclidean distance metric is used.
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept 8.4211e+01 8.4492e+01 8.4726e+01 8.4827e+01 84.8851
## x1 4.9535e-05 5.3409e-05 5.6218e-05 5.7839e-05 0.0001
## x2 -6.0006e-03 -5.6634e-03 -5.2029e-03 -4.8503e-03 -0.0045
## x3 -2.9524e-01 -2.9023e-01 -2.8446e-01 -2.7128e-01 -0.2599
## x4 6.9039e-03 7.1827e-03 7.5466e-03 7.8226e-03 0.0089
## ************************Diagnostic information*************************
## Number of data points: 26
## Effective number of parameters (2trace(S) - trace(S'S)): 5.347672
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 20.65233
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 127.3209
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 115.4373
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 101.1242
## Residual sum of squares: 105.7476
## R-square value: 0.1095733
## Adjusted R-square value: -0.1327242
##
## ***********************************************************************
## Program stops at: 2024-11-28 23:46:20.346394
results4 <-as.data.frame(gwr.model4$SDF)
head(results4)
## Intercept x1 x2 x3 x4 y yhat
## 1 84.63166 5.539258e-05 -0.005075173 -0.2797527 0.007361837 83.090 80.64020
## 2 84.76392 5.727003e-05 -0.005369046 -0.2863261 0.007100894 81.220 80.08621
## 3 84.26299 4.953547e-05 -0.004457572 -0.2610090 0.008122289 79.220 82.05553
## 4 84.88457 6.111864e-05 -0.005895573 -0.2947521 0.007555529 83.140 82.22175
## 5 84.84911 5.790629e-05 -0.005893059 -0.2927155 0.007823111 79.460 81.06788
## 6 84.24865 4.972745e-05 -0.004501894 -0.2606717 0.008301476 80.145 81.27381
## residual CV_Score Stud_residual Intercept_SE x1_SE x2_SE
## 1 2.4498043 0 1.1872840 6.316927 0.0002638607 0.01701724
## 2 1.1337898 0 0.5614789 6.317271 0.0002638522 0.01701742
## 3 -2.8355335 0 -1.4296869 6.318487 0.0002639183 0.01702109
## 4 0.9182457 0 0.4542282 6.318622 0.0002638937 0.01702055
## 5 -1.6078842 0 -0.9847637 6.318530 0.0002639144 0.01702003
## 6 -1.1288070 0 -0.5216459 6.318567 0.0002639155 0.01702037
## x3_SE x4_SE Intercept_TV x1_TV x2_TV x3_TV x4_TV
## 1 0.3121270 0.03782171 13.39760 0.2099311 -0.2982372 -0.8962785 0.1946458
## 2 0.3121413 0.03782374 13.41781 0.2170535 -0.3155030 -0.9172965 0.1877365
## 3 0.3121954 0.03783228 13.33595 0.1876924 -0.2618853 -0.8360439 0.2146920
## 4 0.3122134 0.03782998 13.43403 0.2316033 -0.3463796 -0.9440727 0.1997233
## 5 0.3121870 0.03783331 13.42862 0.2194131 -0.3462426 -0.9376288 0.2067784
## 6 0.3121985 0.03783246 13.33351 0.1884219 -0.2645003 -0.8349549 0.2194273
## Local_R2
## 1 0.1080547
## 2 0.1109331
## 3 0.1008400
## 4 0.1167930
## 5 0.1170374
## 6 0.1012037
bw5 <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar, coords = CoordK[c(1:27),c(1,2)])
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords =
## CoordK[c(1:27), : data is Spatial* object, ignoring coords argument
## Bandwidth: 95.40776 CV score: 168.8706
## Bandwidth: 154.2189 CV score: 175.2807
## Bandwidth: 59.06049 CV score: 174.6529
## Bandwidth: 105.4689 CV score: 169.4827
## Bandwidth: 94.01905 CV score: 168.8408
## Bandwidth: 91.60058 CV score: 168.8292
## Bandwidth: 92.25716 CV score: 168.827
## Bandwidth: 92.28323 CV score: 168.827
## Bandwidth: 92.27427 CV score: 168.827
## Bandwidth: 92.27431 CV score: 168.827
## Bandwidth: 92.27435 CV score: 168.827
## Bandwidth: 92.27431 CV score: 168.827
bw5
## [1] 92.27431
gwr.model5 = gwr(Y0 ~ x1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
bandwidth = bw5,
hatmatrix=TRUE,
se.fit=TRUE) ; gwr.model5
## Warning in gwr(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], bandwidth = bw5, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 92.27431
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 7.8064e+01 8.2843e+01 8.7229e+01 9.1685e+01 9.4016e+01
## x1 -1.5221e-04 -3.5255e-06 9.8135e-05 1.3253e-04 1.5274e-04
## x2 -2.3079e-02 -1.6884e-02 -6.0931e-03 3.2025e-03 1.1024e-02
## x3 -7.4566e-01 -6.3475e-01 -3.8658e-01 -1.1770e-01 1.2204e-01
## x4 -1.8489e-02 -1.0485e-02 -8.1022e-03 -2.6072e-03 1.0435e-02
## Global
## X.Intercept. 84.5552
## x1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 9.538932
## Effective degrees of freedom (residual: 2traceS - traceS'S): 16.46107
## Sigma (residual: 2traceS - traceS'S): 2.008031
## Effective number of parameters (model: traceS): 7.835254
## Effective degrees of freedom (model: traceS): 18.16475
## Sigma (model: traceS): 1.911546
## Sigma (ML): 1.597764
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 126.5742
## AIC (GWR p. 96, eq. 4.22): 105.9875
## Residual sum of squares: 66.37412
## Quasi-global R2: 0.4411099
results5 <-as.data.frame(gwr.model5$SDF)
head(results5)
## sum.w X.Intercept. x1 x2 x3 x4
## 1 17.60304 85.64275 7.266063e-05 -0.003328157 -0.3119534 -0.002536102
## 2 17.88570 88.44540 1.128928e-04 -0.008996573 -0.4549301 -0.007279586
## 3 11.88313 78.06798 -1.522084e-04 0.010960825 0.1220398 0.008799682
## 4 13.68349 93.57458 1.525743e-04 -0.020800124 -0.7063672 -0.014899632
## 5 13.19426 93.66587 1.379342e-04 -0.022160559 -0.7333450 -0.008636579
## 6 12.28392 78.06442 -1.477018e-04 0.011023906 0.1144723 0.010081931
## X.Intercept._se x1_se x2_se x3_se x4_se gwr.e
## 1 5.423000 0.0002300258 0.01490034 0.2678105 0.03260218 2.2346957
## 2 5.532710 0.0002263921 0.01495079 0.2710718 0.03318470 1.0818835
## 3 6.180502 0.0003083922 0.01768941 0.3168566 0.03761820 -0.9876861
## 4 6.340025 0.0002364614 0.01616808 0.3052764 0.03706232 -0.4877132
## 5 6.465809 0.0002379221 0.01638880 0.3055711 0.03836696 -1.3704427
## 6 6.137229 0.0003006047 0.01738543 0.3129749 0.03723194 -0.5994649
## pred pred.se localR2 X.Intercept._se_EDF x1_se_EDF x2_se_EDF
## 1 80.85530 0.7947677 0.3820318 5.696725 0.0002416363 0.01565243
## 2 80.13812 0.8581217 0.4332161 5.811973 0.0002378192 0.01570543
## 3 80.20769 1.0666646 0.2140749 6.492461 0.0003239582 0.01858228
## 4 83.62771 0.9559429 0.5443604 6.660036 0.0002483968 0.01698416
## 5 80.83044 1.4412400 0.5549786 6.792170 0.0002499312 0.01721602
## 6 80.74446 0.6356685 0.2242256 6.447005 0.0003157776 0.01826296
## x3_se_EDF x4_se_EDF pred.se.1
## 1 0.2813282 0.03424777 0.8348834
## 2 0.2847541 0.03485969 0.9014352
## 3 0.3328499 0.03951697 1.1205043
## 4 0.3206852 0.03893304 1.0041939
## 5 0.3209948 0.04030353 1.5139864
## 6 0.3287722 0.03911121 0.6677538
bw6 <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar,
coords = CoordK[c(1:27),c(1,2)], gweight = gwr.bisquare)
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords =
## CoordK[c(1:27), : data is Spatial* object, ignoring coords argument
## Bandwidth: 95.40776 CV score: 180.2693
## Bandwidth: 154.2189 CV score: 184.7979
## Bandwidth: 59.06049 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 117.8716 CV score: 188.9976
## Bandwidth: 81.52434 CV score: 188.7706
## Bandwidth: 103.9882 CV score: 186.3523
## Bandwidth: 90.10477 CV score: 176.5574
## Bandwidth: 86.82734 CV score: 177.4142
## Bandwidth: 89.63268 CV score: 176.4504
## Bandwidth: 89.21726 CV score: 176.4101
## Bandwidth: 89.0938 CV score: 176.4089
## Bandwidth: 89.12484 CV score: 176.4087
## Bandwidth: 89.12581 CV score: 176.4087
## Bandwidth: 89.12572 CV score: 176.4087
## Bandwidth: 89.12576 CV score: 176.4087
## Bandwidth: 89.12568 CV score: 176.4087
## Bandwidth: 89.12572 CV score: 176.4087
bw6
## [1] 89.12572
gwr.model6<-gwr(Y0 ~ x1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
bandwidth = bw6,
gweight = gwr.bisquare,
se.fit=T, hatmatrix=T) ; gwr.model6
## Warning in gwr(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], bandwidth = bw6, gweight = gwr.bisquare, hatmatrix = T,
## se.fit = T)
## Kernel function: gwr.bisquare
## Fixed bandwidth: 89.12572
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 6.2386e+01 8.5749e+01 9.5026e+01 1.0187e+02 1.0632e+02
## x1 -1.0127e-03 2.0830e-04 3.2951e-04 5.5884e-04 1.4769e-03
## x2 -8.4500e-02 -5.0468e-02 -2.3260e-02 -1.2215e-02 6.1570e-02
## x3 -1.4976e+00 -9.9181e-01 -7.4932e-01 8.5500e-02 1.1341e+00
## x4 -1.8666e-01 -6.5184e-02 -3.3328e-02 2.5648e-02 1.1076e-01
## Global
## X.Intercept. 84.5552
## x1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 20.37614
## Effective degrees of freedom (residual: 2traceS - traceS'S): 5.623863
## Sigma (residual: 2traceS - traceS'S): 1.746461
## Effective number of parameters (model: traceS): 17.66862
## Effective degrees of freedom (model: traceS): 8.331377
## Sigma (model: traceS): 1.434888
## Sigma (ML): 0.81225
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 216.2982
## AIC (GWR p. 96, eq. 4.22): 80.64018
## Residual sum of squares: 17.1535
## Quasi-global R2: 0.8555623
results6 <-as.data.frame(gwr.model6$SDF)
results6
## sum.w X.Intercept. x1 x2 x3 x4
## 1 5.949094 97.68681 5.199341e-04 -0.01387242 -1.0002443 -0.03236676
## 2 6.470816 102.61647 6.095857e-04 -0.02359593 -1.1948393 -0.05314593
## 3 4.168484 67.09810 -3.962103e-04 0.01558579 0.8354694 0.03340610
## 4 4.719827 92.96627 2.057689e-04 -0.04179592 -0.7490314 0.02115153
## 5 5.163471 94.73444 3.112166e-04 -0.05510583 -0.8677972 0.02559330
## 6 4.285248 67.53752 -7.013153e-04 0.02314317 0.8248671 0.02995842
## 7 4.606632 105.90521 6.716235e-04 -0.03294798 -1.3370128 -0.07663509
## 8 6.569650 102.80632 6.338338e-04 -0.02292427 -1.1952314 -0.05679384
## 9 4.469891 106.22692 1.465620e-03 -0.06932412 -0.7800593 -0.18423401
## 10 5.086996 91.74038 2.348171e-04 -0.04903682 -0.6885748 0.03306234
## 11 4.287630 62.38596 -1.012668e-03 0.06156978 0.7798006 0.11075803
## 12 3.613063 84.82101 1.122601e-03 -0.08450012 1.1340909 -0.18666374
## 13 4.348704 98.47391 2.360011e-04 -0.03323990 -0.8758367 -0.03428824
## 14 6.555447 102.23706 5.718077e-04 -0.02124618 -1.1183194 -0.06547061
## 15 5.115655 92.47312 2.158980e-04 -0.04556102 -0.7155855 0.02566594
## 16 5.397185 95.31826 3.251248e-04 -0.05597773 -0.8954593 0.02270675
## 17 4.442306 106.04035 1.476854e-03 -0.07059566 -0.7489324 -0.18593386
## 18 4.966038 99.44418 3.508640e-04 -0.05094547 -1.0641936 -0.01019341
## 19 5.073792 67.11569 -7.390496e-04 0.02606396 0.8222382 0.03395609
## 20 5.153757 81.41590 3.839184e-04 -0.02066674 0.1763809 -0.03478512
## 21 5.583159 88.53220 4.498703e-04 -0.02270737 -0.1871412 -0.06432484
## 22 6.377163 100.64682 3.338957e-04 -0.01982826 -0.9665162 -0.06263673
## 23 4.944458 100.78173 4.447228e-05 -0.01166207 -0.7495990 -0.10593708
## 24 3.139827 106.32201 4.520561e-04 -0.06526466 -1.4976279 -0.01299588
## 25 6.095098 94.01097 2.286223e-04 -0.01020872 -0.4954600 -0.07989641
## 26 4.100044 65.46914 -8.539512e-04 0.02930412 0.9226112 0.03801295
## X.Intercept._se x1_se x2_se x3_se x4_se gwr.e
## 1 6.153900 0.0003284358 0.01917437 0.3203081 0.04088139 0.45595345
## 2 6.636673 0.0003005836 0.01784801 0.3453698 0.03719718 0.71494834
## 3 8.440356 0.0007145564 0.03961038 0.5711485 0.04138352 -0.02890283
## 4 9.063194 0.0002288521 0.01956415 0.3728415 0.06153037 -0.92128964
## 5 9.180487 0.0002655403 0.02232960 0.4242468 0.05345068 0.10658034
## 6 7.463207 0.0005873769 0.03332000 0.4704016 0.04569217 -0.98705048
## 7 7.540380 0.0003477401 0.02207048 0.3759970 0.04644213 -1.67657938
## 8 7.021027 0.0003328734 0.01881807 0.3623992 0.04008081 0.21481565
## 9 20.383772 0.0007575030 0.03789628 1.1341898 0.10025772 0.08577090
## 10 9.316596 0.0002457771 0.02047170 0.4400542 0.05453730 0.14222500
## 11 7.300013 0.0004436540 0.02295552 0.4078406 0.04638133 -0.32494804
## 12 10.387730 0.0006543649 0.04391819 0.8755532 0.07592018 0.14556422
## 13 7.850050 0.0002213641 0.01834086 0.3399128 0.04706456 0.41584745
## 14 7.105043 0.0003321778 0.01941967 0.3522365 0.04224892 0.48509815
## 15 8.634476 0.0002344485 0.01960322 0.3905381 0.05478593 0.34554149
## 16 8.506282 0.0002521863 0.02135292 0.3867945 0.05057492 0.46759797
## 17 20.468402 0.0007611624 0.03817679 1.1416923 0.10088653 0.34577537
## 18 8.252017 0.0003481871 0.02815831 0.3936712 0.04376632 -0.25872959
## 19 7.100636 0.0005208575 0.02545144 0.4407431 0.04405998 0.65243740
## 20 10.988296 0.0005899686 0.03209844 0.6875416 0.05001804 2.38355755
## 21 8.355464 0.0004733351 0.02347956 0.4955265 0.04052701 -1.87306253
## 22 6.075333 0.0002165391 0.01422600 0.2855350 0.03510698 -0.10944578
## 23 6.991961 0.0003296897 0.01734277 0.3299176 0.04125638 0.57101117
## 24 19.722087 0.0004912461 0.04569061 0.9731700 0.10674485 -0.21365102
## 25 5.733123 0.0003042213 0.01691750 0.2940518 0.03393143 0.47877784
## 26 7.506626 0.0005923430 0.03223079 0.4553513 0.05095034 -0.65322343
## pred pred.se localR2 X.Intercept._se_EDF x1_se_EDF x2_se_EDF
## 1 82.63405 1.2285540 0.8999262 7.490167 0.0003997528 0.02333792
## 2 80.50505 0.8735137 0.8906165 8.077769 0.0003658527 0.02172355
## 3 79.24890 1.2734963 0.5770913 10.273106 0.0008697161 0.04821143
## 4 84.06129 1.0919966 0.9393527 11.031188 0.0002785453 0.02381233
## 5 79.35342 1.3734723 0.9639721 11.173950 0.0003232000 0.02717828
## 6 81.13205 1.0571816 0.6413588 9.083777 0.0007149207 0.04055515
## 7 79.49658 0.8913851 0.8453972 9.177709 0.0004232488 0.02686289
## 8 83.63518 1.2755990 0.8957206 8.545583 0.0004051540 0.02290425
## 9 80.63423 1.3168422 0.2903513 24.809931 0.0009219882 0.04612513
## 10 84.34777 1.0628915 0.9588111 11.339614 0.0002991455 0.02491696
## 11 80.09495 1.2254962 0.9362090 8.885147 0.0005399896 0.02794011
## 12 80.08444 0.9636901 0.5086527 12.643335 0.0007964546 0.05345465
## 13 80.53415 1.1480248 0.9296095 9.554620 0.0002694314 0.02232342
## 14 83.74490 1.3038214 0.8923839 8.647841 0.0004043073 0.02363648
## 15 85.00446 0.9811267 0.9472961 10.509377 0.0002853569 0.02385989
## 16 82.81240 0.7741018 0.9596191 10.353348 0.0003069463 0.02598952
## 17 80.79422 1.2128204 0.2900807 24.912938 0.0009264422 0.04646655
## 18 81.56873 1.1330836 0.8850742 10.043871 0.0004237929 0.03427265
## 19 80.31756 1.1256649 0.7127814 8.642478 0.0006339572 0.03097800
## 20 80.73644 0.8350579 0.3680830 13.374309 0.0007180752 0.03906833
## 21 79.96306 0.8920217 0.4908546 10.169780 0.0005761157 0.02857794
## 22 79.28945 1.0601545 0.9174277 7.394538 0.0002635587 0.01731505
## 23 80.84899 0.7126806 0.8863097 8.510205 0.0004012790 0.02110861
## 24 79.50365 1.1053694 0.8443847 24.004567 0.0005979159 0.05561193
## 25 80.92122 0.5319267 0.7867199 6.978021 0.0003702803 0.02059099
## 26 77.50322 1.3286808 0.8026318 9.136624 0.0007209651 0.03922942
## x3_se_EDF x4_se_EDF pred.se.1
## 1 0.3898602 0.04975843 1.4953238
## 2 0.4203638 0.04527422 1.0631896
## 3 0.6951684 0.05036960 1.5500250
## 4 0.4538008 0.07489115 1.3291141
## 5 0.5163683 0.06505703 1.6717099
## 6 0.5725452 0.05561383 1.2867394
## 7 0.4576415 0.05652664 1.0849416
## 8 0.4410910 0.04878401 1.5525842
## 9 1.3804693 0.12202781 1.6027830
## 10 0.5356081 0.06637960 1.2936891
## 11 0.4963996 0.05645263 1.4916020
## 12 1.0656719 0.09240558 1.1729471
## 13 0.4137219 0.05728422 1.3973084
## 14 0.4287216 0.05142291 1.5869348
## 15 0.4753401 0.06668222 1.1941699
## 16 0.4707836 0.06155683 0.9421913
## 17 1.3896009 0.12279317 1.4761738
## 18 0.4791535 0.05326979 1.3791228
## 19 0.5364466 0.05362722 1.3700933
## 20 0.8368353 0.06087902 1.0163834
## 21 0.6031258 0.04932710 1.0857165
## 22 0.3475364 0.04273016 1.2903578
## 23 0.4015564 0.05021485 0.8674331
## 24 1.1844854 0.12992357 1.3453908
## 25 0.3579026 0.04129934 0.6474299
## 26 0.5542269 0.06201377 1.6171923
bw7 <- gwr.sel(Y0 ~ x1+x2+x3+x4, data = Jabar,
coords = CoordK[c(1:27),c(1,2)], gweight = gwr.tricube)
## Warning in gwr.sel(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords =
## CoordK[c(1:27), : data is Spatial* object, ignoring coords argument
## Bandwidth: 95.40776 CV score: 177.3272
## Bandwidth: 154.2189 CV score: 189.771
## Bandwidth: 59.06049 CV score: NA
## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value
## Bandwidth: 117.8716 CV score: 190.8357
## Bandwidth: 81.52434 CV score: 206.6929
## Bandwidth: 103.9882 CV score: 184.1265
## Bandwidth: 90.10477 CV score: 177.504
## Bandwidth: 93.03645 CV score: 176.5313
## Bandwidth: 92.88867 CV score: 176.5218
## Bandwidth: 92.66462 CV score: 176.5179
## Bandwidth: 92.70736 CV score: 176.5176
## Bandwidth: 92.70852 CV score: 176.5176
## Bandwidth: 92.70832 CV score: 176.5176
## Bandwidth: 92.70836 CV score: 176.5176
## Bandwidth: 92.70827 CV score: 176.5176
## Bandwidth: 92.70832 CV score: 176.5176
bw7
## [1] 92.70832
gwr.model7<-gwr(Y0 ~ x1+x2+x3+x4,
data = Jabar,
coords=CoordK[c(1:27),c(1,2)],
bandwidth = bw7,
gweight = gwr.tricube,
se.fit=T, hatmatrix=T) ; gwr.model7
## Warning in gwr(Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27), :
## data is Spatial* object, ignoring coords argument
## Call:
## gwr(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, coords = CoordK[c(1:27),
## c(1, 2)], bandwidth = bw7, gweight = gwr.tricube, hatmatrix = T,
## se.fit = T)
## Kernel function: gwr.tricube
## Fixed bandwidth: 92.70832
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max.
## X.Intercept. 6.2425e+01 8.5532e+01 9.4810e+01 1.0146e+02 1.0597e+02
## x1 -1.0195e-03 2.0219e-04 3.3502e-04 5.4840e-04 1.4366e-03
## x2 -8.3936e-02 -5.0834e-02 -2.3146e-02 -1.2789e-02 6.1237e-02
## x3 -1.4799e+00 -9.7351e-01 -7.5357e-01 1.8249e-01 1.1352e+00
## x4 -1.8504e-01 -6.2455e-02 -2.6611e-02 2.5868e-02 1.0948e-01
## Global
## X.Intercept. 84.5552
## x1 0.0001
## x2 -0.0052
## x3 -0.2773
## x4 0.0081
## Number of data points: 26
## Effective number of parameters (residual: 2traceS - traceS'S): 19.06791
## Effective degrees of freedom (residual: 2traceS - traceS'S): 6.932087
## Sigma (residual: 2traceS - traceS'S): 1.761698
## Effective number of parameters (model: traceS): 16.65357
## Effective degrees of freedom (model: traceS): 9.34643
## Sigma (model: traceS): 1.517192
## Sigma (ML): 0.9096549
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 193.8176
## AIC (GWR p. 96, eq. 4.22): 85.51449
## Residual sum of squares: 21.51427
## Quasi-global R2: 0.8188433
results7 <-as.data.frame(gwr.model7$SDF)
bw8 <- bw.gwr(Y0 ~ x1+x2+x3+x4,
data=Jabar,
approach="CV",
kernel="exponential",
adaptive=F)
## Fixed bandwidth: 1.23965 CV score: 170.1828
## Fixed bandwidth: 0.766299 CV score: 165.3876
## Fixed bandwidth: 0.473752 CV score: 163.7527
## Fixed bandwidth: 0.2929481 CV score: 165.6208
## Fixed bandwidth: 0.585495 CV score: 163.9286
## Fixed bandwidth: 0.4046911 CV score: 164.1391
## Fixed bandwidth: 0.5164341 CV score: 163.719
## Fixed bandwidth: 0.542813 CV score: 163.7647
## Fixed bandwidth: 0.500131 CV score: 163.7152
## Fixed bandwidth: 0.4900551 CV score: 163.7229
## Fixed bandwidth: 0.5063582 CV score: 163.7143
## Fixed bandwidth: 0.5102068 CV score: 163.7152
## Fixed bandwidth: 0.5039796 CV score: 163.7143
bw8
## [1] 0.5039796
gwr.model8 <- gwr.basic(Y0 ~ x1+x2+x3+x4,
data=Jabar,
bw=bw8,
kernel="exponential") ; gwr.model8
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2024-11-28 23:46:20.806801
## Call:
## gwr.basic(formula = Y0 ~ x1 + x2 + x3 + x4, data = Jabar, bw = bw8,
## kernel = "exponential")
##
## Dependent (y) variable: Y0
## Independent variables: x1 x2 x3 x4
## Number of data points: 26
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6270 -1.3180 0.0582 1.6168 3.3910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.456e+01 6.359e+00 13.296 1.07e-11 ***
## x1 5.465e-05 2.656e-04 0.206 0.839
## x2 -5.238e-03 1.713e-02 -0.306 0.763
## x3 -2.773e-01 3.142e-01 -0.883 0.387
## x4 8.083e-03 3.807e-02 0.212 0.834
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 2.278 on 21 degrees of freedom
## Multiple R-squared: 0.08228
## Adjusted R-squared: -0.09253
## F-statistic: 0.4707 on 4 and 21 DF, p-value: 0.7566
## ***Extra Diagnostic information
## Residual sum of squares: 108.9892
## Sigma(hat): 2.131013
## AIC: 123.0468
## AICc: 127.4678
## BIC: 124.1439
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: exponential
## Fixed bandwidth: 0.5039796
## Regression points: the same locations as observations are used.
## Distance metric: Euclidean distance metric is used.
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept 7.0613e+01 8.4409e+01 9.1808e+01 9.3845e+01 94.6141
## x1 -3.6010e-04 7.6101e-05 1.0899e-04 1.6736e-04 0.0003
## x2 -3.2414e-02 -2.0586e-02 -7.0291e-03 -3.7787e-04 0.0184
## x3 -7.7949e-01 -7.2635e-01 -5.5711e-01 -1.5344e-01 0.4722
## x4 -3.6922e-02 -2.4616e-02 -1.6094e-02 2.2756e-03 0.0460
## ************************Diagnostic information*************************
## Number of data points: 26
## Effective number of parameters (2trace(S) - trace(S'S)): 17.03576
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 8.964241
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 146.2305
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 92.23924
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 95.78053
## Residual sum of squares: 31.96731
## R-square value: 0.7308256
## Adjusted R-square value: 0.1550531
##
## ***********************************************************************
## Program stops at: 2024-11-28 23:46:20.812178
results8 <-as.data.frame(gwr.model8$SDF)
head(results8)
## Intercept x1 x2 x3 x4 y yhat
## 1 88.23801 0.0001407625 -0.0002064902 -0.4374236 -0.016669285 83.090 81.56384
## 2 93.10286 0.0001905282 -0.0114147488 -0.6437838 -0.032230233 81.220 80.45520
## 3 75.15907 -0.0002699751 0.0168253488 0.2925635 0.014754251 79.220 79.70282
## 4 93.65257 0.0001949110 -0.0285654047 -0.7170841 -0.006339371 83.140 83.78157
## 5 94.09813 0.0001016543 -0.0287306020 -0.7727522 0.002642173 79.460 80.01653
## 6 74.10167 -0.0002734831 0.0162184131 0.3340180 0.021628146 80.145 80.71545
## residual CV_Score Stud_residual Intercept_SE x1_SE x2_SE
## 1 1.5261562 0 1.3696676 5.723367 0.0002603256 0.01670844
## 2 0.7648013 0 0.6330562 6.162243 0.0002413674 0.01676064
## 3 -0.4828211 0 -0.6057993 7.132206 0.0004037802 0.02255295
## 4 -0.6415683 0 -0.6302865 7.366547 0.0002545922 0.01852380
## 5 -0.5565253 0 -1.0018514 7.526387 0.0002524468 0.01730035
## 6 -0.5704542 0 -0.4733771 7.237053 0.0004026689 0.02163847
## x3_SE x4_SE Intercept_TV x1_TV x2_TV x3_TV
## 1 0.2790713 0.03533118 15.41715 0.5407172 -0.01235844 -1.5674258
## 2 0.3009940 0.03670549 15.10860 0.7893703 -0.68104481 -2.1388591
## 3 0.3914284 0.04216437 10.53798 -0.6686190 0.74603762 0.7474255
## 4 0.3625146 0.04056815 12.71323 0.7655809 -1.54209220 -1.9780833
## 5 0.3405634 0.04540292 12.50243 0.4026762 -1.66069500 -2.2690406
## 6 0.3941056 0.04291225 10.23920 -0.6791760 0.74951764 0.8475343
## x4_TV Local_R2
## 1 -0.47180102 0.6784367
## 2 -0.87807664 0.7129397
## 3 0.34992229 0.6520996
## 4 -0.15626473 0.8600502
## 5 0.05819389 0.8560341
## 6 0.50400860 0.6624303
Perbandingan_AIC_GWR = data.frame("Model"=c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6", "Model 7", "Model 8"),
"AIC"=c(gwr.model1$results$AICh,gwr.model2$results$AICh,gwr.model3$results$AICh,gwr.model4$GW.diagnostic$AIC,gwr.model5$results$AICh,gwr.model6$results$AICh,gwr.model7$results$AICh,gwr.model8$GW.diagnostic$AIC))
Perbandingan_AIC_GWR
## Model AIC
## 1 Model 1 105.79763
## 2 Model 2 106.34468
## 3 Model 3 106.80369
## 4 Model 4 115.43728
## 5 Model 5 105.98753
## 6 Model 6 80.64018
## 7 Model 7 85.51449
## 8 Model 8 92.23924
resid_gwr6 <- results6$gwr.e
Jabar$resid_gwr6 <- resid_gwr6
#spplot(Jabar,zcol="resid_gwr6",main="Peta Sebaran Residual Model Fixed Bisquare")
tm_shape(Jabar) +
tm_polygons("resid_gwr6", title = "Residual GWR") +
tm_layout(title = "Peta Sebaran Residual Model Fixed Bisquare")
## Variable(s) "resid_gwr6" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
#Residuals diagnostics
shapiro.test(resid_gwr6) #Normality
##
## Shapiro-Wilk normality test
##
## data: resid_gwr6
## W = 0.90016, p-value = 0.01578
library(lmtest)
bptest(gwr.model6$lm, weights = gwr.model6$gweight) #Heteroskedasticity
##
## studentized Breusch-Pagan test
##
## data: gwr.model6$lm
## BP = 8.2735, df = 4, p-value = 0.08206
gwr.morantest(gwr.model6, W_optimum, zero.policy = TRUE) #Non-autocorrelation
##
## Leung et al. 2000 three moment approximation for Moran's I
##
## data: GWR residuals
## statistic = 9.5755, df = 12.249, p-value = 0.3272
## sample estimates:
## I
## -0.1594028
#Comparing OLS and GWR models under an inferential framework
BFC02.gwr.test(gwr.model6)
##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: gwr.model6
## F = 6.3538, df1 = 21.0000, df2 = 5.6239, p-value = 0.01782
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 108.9892 17.1535
#Parameter Significance Test
LMZ.F3GWR.test(gwr.model6)
##
## Leung et al. (2000) F(3) test
##
## F statistic Numerator d.f. Denominator d.f. Pr(>)
## (Intercept) 2.5635e+08 7.7793e+00 9.8571 < 2.2e-16 ***
## x1 3.0135e+08 1.1823e+01 9.8571 < 2.2e-16 ***
## x2 2.0549e+08 8.4209e+00 9.8571 < 2.2e-16 ***
## x3 4.1267e+08 7.4249e+00 9.8571 < 2.2e-16 ***
## x4 1.6669e+08 8.3223e+00 9.8571 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Comparing AIC
data.frame("MODEL" = c("GWR","OLS"),
"AIC" = c(gwr.model6$results$AICh,AIC(jabarlm)),
"R2" = c(0.8545739,0.07381))%>% arrange(AIC)
## MODEL AIC R2
## 1 GWR 80.64018 0.8545739
## 2 OLS 123.04675 0.0738100
#Examine the correlation
corr <- round(cor(as.data.frame(results4[,2:6]), use ="complete.obs"),2)
Koefisien pada tiap lokasi
data_koefisien <- data.frame(
Nama_Daerah = data$Kabupaten.Kota,
Intersept = results6$X.Intercept.,
Koefisien_X1 = results6$x1, # Koefisien untuk variabel X1
Koefisien_X2 = results6$x2, # Koefisien untuk variabel x2
Koefisien_X3 = results6$x3, # Koefisien untuk variabel x3
Koefisien_X4 = results6$x4 # Koefisien untuk variabel x4
)
data_koefisien
## Nama_Daerah Intersept Koefisien_X1 Koefisien_X2 Koefisien_X3
## 1 Bandung 97.68681 5.199341e-04 -0.01387242 -1.0002443
## 2 Bandung Barat 102.61647 6.095857e-04 -0.02359593 -1.1948393
## 3 Banjar 67.09810 -3.962103e-04 0.01558579 0.8354694
## 4 Bekasi 92.96627 2.057689e-04 -0.04179592 -0.7490314
## 5 Bogor 94.73444 3.112166e-04 -0.05510583 -0.8677972
## 6 Ciamis 67.53752 -7.013153e-04 0.02314317 0.8248671
## 7 Cianjur 105.90521 6.716235e-04 -0.03294798 -1.3370128
## 8 Cimahi 102.80632 6.338338e-04 -0.02292427 -1.1952314
## 9 Cirebon 106.22692 1.465620e-03 -0.06932412 -0.7800593
## 10 Depok 91.74038 2.348171e-04 -0.04903682 -0.6885748
## 11 Garut 62.38596 -1.012668e-03 0.06156978 0.7798006
## 12 Indramayu 84.82101 1.122601e-03 -0.08450012 1.1340909
## 13 Karawang 98.47391 2.360011e-04 -0.03323990 -0.8758367
## 14 Kota Bandung 102.23706 5.718077e-04 -0.02124618 -1.1183194
## 15 Kota Bekasi 92.47312 2.158980e-04 -0.04556102 -0.7155855
## 16 Kota Bogor 95.31826 3.251248e-04 -0.05597773 -0.8954593
## 17 Kota Cirebon 106.04035 1.476854e-03 -0.07059566 -0.7489324
## 18 Kota Sukabumi 99.44418 3.508640e-04 -0.05094547 -1.0641936
## 19 Kota Tasikmalaya 67.11569 -7.390496e-04 0.02606396 0.8222382
## 20 Kuningan 81.41590 3.839184e-04 -0.02066674 0.1763809
## 21 Majalengka 88.53220 4.498703e-04 -0.02270737 -0.1871412
## 22 Purwakarta 100.64682 3.338957e-04 -0.01982826 -0.9665162
## 23 Subang 100.78173 4.447228e-05 -0.01166207 -0.7495990
## 24 Sukabumi 106.32201 4.520561e-04 -0.06526466 -1.4976279
## 25 Sumedang 94.01097 2.286223e-04 -0.01020872 -0.4954600
## 26 Tasikmalaya 65.46914 -8.539512e-04 0.02930412 0.9226112
## Koefisien_X4
## 1 -0.03236676
## 2 -0.05314593
## 3 0.03340610
## 4 0.02115153
## 5 0.02559330
## 6 0.02995842
## 7 -0.07663509
## 8 -0.05679384
## 9 -0.18423401
## 10 0.03306234
## 11 0.11075803
## 12 -0.18666374
## 13 -0.03428824
## 14 -0.06547061
## 15 0.02566594
## 16 0.02270675
## 17 -0.18593386
## 18 -0.01019341
## 19 0.03395609
## 20 -0.03478512
## 21 -0.06432484
## 22 -0.06263673
## 23 -0.10593708
## 24 -0.01299588
## 25 -0.07989641
## 26 0.03801295
#Map the results
gwr.map <- cbind(Jabar, as.matrix(results6))
qtm(gwr.map, fill = "localR2")
gwr.map2 <- st_as_sf(gwr.map)
#Local significance test
t1 = gwr.model6$SDF$x1 / gwr.model6$SDF$x1_se
t2 = gwr.model6$SDF$x2 / gwr.model6$SDF$x2_se
t3 = gwr.model6$SDF$x3 / gwr.model6$SDF$x3_se
t4 = gwr.model6$SDF$x4 / gwr.model6$SDF$x4_se
pvalue1 = 2*pt(q=t1, df=25, lower.tail=TRUE)
pvalue2 = 2*pt(q=t2, df=25, lower.tail=TRUE)
pvalue3 = 2*pt(q=t3, df=25, lower.tail=TRUE)
pvalue4 = 2*pt(q=t4, df=25, lower.tail=TRUE)
data.frame(Kabupaten.Kota, t1, pvalue1) %>% mutate(ket=ifelse(pvalue1<0.05,"Reject", "Not Reject"))
## Kabupaten.Kota t1 pvalue1 ket
## 1 Bandung 1.5830615 1.87402255 Not Reject
## 2 Bandung Barat 2.0280076 1.94665748 Not Reject
## 3 Banjar -0.5544844 0.58417428 Not Reject
## 4 Bekasi 0.8991349 1.62284009 Not Reject
## 5 Bogor 1.1720129 1.74775977 Not Reject
## 6 Ciamis -1.1939783 0.24369307 Not Reject
## 7 Cianjur 1.9313952 1.93515344 Not Reject
## 8 Cimahi 1.9041287 1.93153616 Not Reject
## 9 Cirebon 1.9348044 1.93559376 Not Reject
## 10 Depok 0.9554066 1.65147993 Not Reject
## 11 Garut -2.2825618 0.03122939 Reject
## 12 Indramayu 1.7155582 1.90138697 Not Reject
## 13 Karawang 1.0661218 1.70344312 Not Reject
## 14 Kota Bandung 1.7213904 1.90246605 Not Reject
## 15 Kota Bekasi 0.9208761 1.63408345 Not Reject
## 16 Kota Bogor 1.2892247 1.79087865 Not Reject
## 17 Kota Cirebon 1.9402609 1.93629303 Not Reject
## 18 Kota Sukabumi 1.0076882 1.67674152 Not Reject
## 19 Kota Tasikmalaya -1.4189092 0.16827337 Not Reject
## 20 Kuningan 0.6507437 1.47885148 Not Reject
## 21 Majalengka 0.9504266 1.64900593 Not Reject
## 22 Purwakarta 1.5419645 1.86435357 Not Reject
## 23 Subang 0.1348913 1.10622201 Not Reject
## 24 Sukabumi 0.9202232 1.63374908 Not Reject
## 25 Sumedang 0.7515000 1.54063160 Not Reject
## 26 Tasikmalaya -1.4416498 0.16181157 Not Reject
data.frame(Kabupaten.Kota, t2, pvalue2) %>% mutate(ket=ifelse(pvalue2<0.05,"Reject", "Not Reject"))
## Kabupaten.Kota t2 pvalue2 ket
## 1 Bandung -0.7234876 0.47609587 Not Reject
## 2 Bandung Barat -1.3220483 0.19812163 Not Reject
## 3 Banjar 0.3934775 1.30269742 Not Reject
## 4 Bekasi -2.1363528 0.04262467 Reject
## 5 Bogor -2.4678375 0.02078753 Reject
## 6 Ciamis 0.6945729 1.50627201 Not Reject
## 7 Cianjur -1.4928528 0.14798949 Not Reject
## 8 Cimahi -1.2182053 0.23451828 Not Reject
## 9 Cirebon -1.8293120 0.07930561 Not Reject
## 10 Depok -2.3953462 0.02441566 Reject
## 11 Garut 2.6821339 1.98722485 Not Reject
## 12 Indramayu -1.9240345 0.06580620 Not Reject
## 13 Karawang -1.8123415 0.08196160 Not Reject
## 14 Kota Bandung -1.0940546 0.28436074 Not Reject
## 15 Kota Bekasi -2.3241595 0.02853650 Reject
## 16 Kota Bogor -2.6215492 0.01468421 Reject
## 17 Kota Cirebon -1.8491775 0.07629115 Not Reject
## 18 Kota Sukabumi -1.8092516 0.08245334 Not Reject
## 19 Kota Tasikmalaya 1.0240662 1.68438832 Not Reject
## 20 Kuningan -0.6438550 0.52553273 Not Reject
## 21 Majalengka -0.9671125 0.34275121 Not Reject
## 22 Purwakarta -1.3938047 0.17564372 Not Reject
## 23 Subang -0.6724454 0.50746823 Not Reject
## 24 Sukabumi -1.4284041 0.16555078 Not Reject
## 25 Sumedang -0.6034414 0.55165138 Not Reject
## 26 Tasikmalaya 0.9091964 1.62807118 Not Reject
data.frame(Kabupaten.Kota, t3, pvalue3) %>% mutate(ket=ifelse(pvalue3<0.05,"Reject", "Not Reject"))
## Kabupaten.Kota t3 pvalue3 ket
## 1 Bandung -3.1227568 0.004487355 Reject
## 2 Bandung Barat -3.4595942 0.001953514 Reject
## 3 Banjar 1.4627884 1.844015378 Not Reject
## 4 Bekasi -2.0089807 0.055453962 Not Reject
## 5 Bogor -2.0455006 0.051464442 Not Reject
## 6 Ciamis 1.7535379 1.908235835 Not Reject
## 7 Cianjur -3.5559134 0.001534120 Reject
## 8 Cimahi -3.2981074 0.002918841 Reject
## 9 Cirebon -0.6877678 0.497931049 Not Reject
## 10 Depok -1.5647500 0.130213022 Not Reject
## 11 Garut 1.9120231 1.932601147 Not Reject
## 12 Indramayu 1.2952850 1.792944178 Not Reject
## 13 Karawang -2.5766516 0.016267549 Reject
## 14 Kota Bandung -3.1749105 0.003951341 Reject
## 15 Kota Bekasi -1.8323065 0.078844735 Not Reject
## 16 Kota Bogor -2.3150776 0.029105673 Reject
## 17 Kota Cirebon -0.6559845 0.517826467 Not Reject
## 18 Kota Sukabumi -2.7032547 0.012166213 Reject
## 19 Kota Tasikmalaya 1.8655726 1.926121562 Not Reject
## 20 Kuningan 0.2565385 1.200365559 Not Reject
## 21 Majalengka -0.3776614 0.708869547 Not Reject
## 22 Purwakarta -3.3849307 0.002353466 Reject
## 23 Subang -2.2720793 0.031943535 Reject
## 24 Sukabumi -1.5389171 0.136387004 Not Reject
## 25 Sumedang -1.6849412 0.104444909 Not Reject
## 26 Tasikmalaya 2.0261523 1.946454776 Not Reject
data.frame(Kabupaten.Kota, t4, pvalue4) %>% mutate(ket=ifelse(pvalue4<0.05,"Reject", "Not Reject"))
## Kabupaten.Kota t4 pvalue4 ket
## 1 Bandung -0.7917237 0.43596991 Not Reject
## 2 Bandung Barat -1.4287624 0.16544873 Not Reject
## 3 Banjar 0.8072320 1.57285340 Not Reject
## 4 Bekasi 0.3437576 1.26609699 Not Reject
## 5 Bogor 0.4788209 1.36377256 Not Reject
## 6 Ciamis 0.6556577 1.48196674 Not Reject
## 7 Cianjur -1.6501200 0.11142798 Not Reject
## 8 Cimahi -1.4169832 0.16882995 Not Reject
## 9 Cirebon -1.8376042 0.07803505 Not Reject
## 10 Depok 0.6062335 1.45017466 Not Reject
## 11 Garut 2.3879876 1.97518519 Not Reject
## 12 Indramayu -2.4586842 0.02121646 Reject
## 13 Karawang -0.7285364 0.47305511 Not Reject
## 14 Kota Bandung -1.5496396 0.13379583 Not Reject
## 15 Kota Bekasi 0.4684768 1.35649725 Not Reject
## 16 Kota Bogor 0.4489726 1.34268070 Not Reject
## 17 Kota Cirebon -1.8429998 0.07721778 Not Reject
## 18 Kota Sukabumi -0.2329053 0.81773124 Not Reject
## 19 Kota Tasikmalaya 0.7706788 1.55188004 Not Reject
## 20 Kuningan -0.6954514 0.49318689 Not Reject
## 21 Majalengka -1.5872091 0.12503406 Not Reject
## 22 Purwakarta -1.7841675 0.08653997 Not Reject
## 23 Subang -2.5677743 0.01659888 Reject
## 24 Sukabumi -0.1217471 0.90407264 Not Reject
## 25 Sumedang -2.3546434 0.02669960 Reject
## 26 Tasikmalaya 0.7460784 1.53742160 Not Reject
melihat persebaran nilai koefisien di tiap kabupaten/kota
library(RColorBrewer)
spplot(gwr.model6$SDF, "x1", main = "Koefisien GWR untuk x1", col.regions = colorRampPalette(brewer.pal(9, "YlGnBu"))(100))
spplot(gwr.model6$SDF, "x2", main = "Koefisien GWR untuk x2", col.regions = colorRampPalette(brewer.pal(9, "YlGnBu"))(100))
spplot(gwr.model6$SDF, "x3", main = "Koefisien GWR untuk x3", col.regions = colorRampPalette(brewer.pal(9, "YlGnBu"))(100))
spplot(gwr.model6$SDF, "x4", main = "Koefisien GWR untuk x4", col.regions = colorRampPalette(brewer.pal(9, "YlGnBu"))(100))