library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(maps)
## Warning: package 'maps' was built under R version 4.3.3
library(raster)
## Warning: package 'raster' was built under R version 4.3.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.3
library(sp)
library(spdep)
## Warning: package 'spdep' was built under R version 4.3.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.3.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(gstat)
## Warning: package 'gstat' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(sf)
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.3.3
## 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
data<-read.csv("D:/Matkul Sem 5/Spasial/uts spasial.csv", sep=";", dec=',')
Indo<-readRDS('D:/Matkul Sem 5/Spasial/gadm36_IDN_2_sp.rds')
Jatim <- Indo[Indo$NAME_1 == "Jawa Timur",]
Jatim$id <- c(1:38)
file_Jatim <- st_as_sf(Jatim)
data$id <- seq_len(nrow(data))
file_Jatim$id <- as.integer(file_Jatim$id)
data$id <- as.integer(data$id)
Jatim_merged <- file_Jatim %>% left_join(data, by="id")
names(Jatim_merged) <- abbreviate(names(Jatim_merged), minlength = 10)
library(psych)
## Warning: package 'psych' was built under R version 4.3.3
##
## Attaching package: 'psych'
## The following object is masked from 'package:raster':
##
## distance
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(spdep)
library(raster)
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.2
## corrplot 0.92 loaded
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
##
## Attaching package: 'DescTools'
## The following objects are masked from 'package:psych':
##
## AUC, ICC, SD
library(nortest)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
##
## Recode
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
library(spatialreg)
library(sf)
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ✔ readr 2.1.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%() masks ggplot2::%+%()
## ✖ psych::alpha() masks ggplot2::alpha()
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::map() masks maps::map()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ car::recode() masks dplyr::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
library(spgwr)
## Warning: package 'spgwr' was built under R version 4.3.3
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
Mapping
Jatim_centroid <- st_centroid(Jatim_merged)
## Warning: st_centroid assumes attributes are constant over geometries
ggplot() +
geom_sf(data = Jatim_merged, aes(fill = Jmlh.K.TBC), color = "black") +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
# Menambahkan nama kabupaten
geom_sf_text(data = Jatim_centroid, aes(label = Kabptn..Kt), size = 1.5, color = "black") +
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 = "Jumlah Kasus TBC")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Statistik Deskriptif
# Check summary of the data
summary(data)
## Kabupaten..Kota Jumlah.Kasus.TBC
## Length:38 Min. : 157.0
## Class :character 1st Qu.: 574.2
## Mode :character Median : 912.0
## Mean :1097.2
## 3rd Qu.:1404.2
## Max. :4101.0
## Persentase.TTU.yang.memenuhi.Syarat.Kesehatan Jumlah.Puskesmas
## Min. :19.00 Min. : 3.00
## 1st Qu.:62.00 1st Qu.:20.25
## Median :69.00 Median :25.00
## Mean :67.89 Mean :25.55
## 3rd Qu.:77.75 3rd Qu.:33.00
## Max. :95.00 Max. :63.00
## Tingkat.Pengangguran.Terbuka
## Min. : 2.000
## 1st Qu.: 4.000
## Median : 5.000
## Mean : 5.579
## 3rd Qu.: 6.750
## Max. :11.000
describe(data)
## vars n mean sd median
## Kabupaten..Kota* 1 38 19.50 11.11 19.5
## Jumlah.Kasus.TBC 2 38 1097.18 786.50 912.0
## Persentase.TTU.yang.memenuhi.Syarat.Kesehatan 3 38 67.89 13.99 69.0
## Jumlah.Puskesmas 4 38 25.55 12.86 25.0
## Tingkat.Pengangguran.Terbuka 5 38 5.58 2.07 5.0
## trimmed mad min max range
## Kabupaten..Kota* 19.50 14.08 1 38 37
## Jumlah.Kasus.TBC 989.12 630.85 157 4101 3944
## Persentase.TTU.yang.memenuhi.Syarat.Kesehatan 68.94 11.86 19 95 76
## Jumlah.Puskesmas 24.97 11.12 3 63 60
## Tingkat.Pengangguran.Terbuka 5.41 1.48 2 11 9
## skew kurtosis se
## Kabupaten..Kota* 0.00 -1.30 1.80
## Jumlah.Kasus.TBC 1.71 3.73 127.59
## Persentase.TTU.yang.memenuhi.Syarat.Kesehatan -1.07 2.15 2.27
## Jumlah.Puskesmas 0.35 0.47 2.09
## Tingkat.Pengangguran.Terbuka 0.80 0.17 0.34
Metode Ordinary Least Square
#Define datasets
Y <- data$Jumlah.Kasus.TBC
X1 <- data$Persentase.TTU.yang.memenuhi.Syarat.Kesehatan
X2 <- data$Jumlah.Puskesmas
X3 <- data$Tingkat.Pengangguran.Terbuka
# Linear Regression (OLS)
model_ols <- lm(Y ~ X1 + X2 + X3, data=data)
summary(model_ols)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -527.05 -159.43 7.16 190.16 477.86
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -551.095 295.170 -1.867 0.0705 .
## X1 -8.707 3.888 -2.240 0.0318 *
## X2 48.549 3.880 12.514 2.79e-14 ***
## X3 179.043 24.322 7.361 1.56e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 278.9 on 34 degrees of freedom
## Multiple R-squared: 0.8844, Adjusted R-squared: 0.8742
## F-statistic: 86.72 on 3 and 34 DF, p-value: 5.263e-16
Uji Normalitas
#Normality test for residuals
residuals <- residuals(model_ols)
shapiro.test(residuals)
##
## Shapiro-Wilk normality test
##
## data: residuals
## W = 0.97999, p-value = 0.7173
qqnorm(model_ols$residuals)
qqline(model_ols$residuals, col="red")
Uji Multikolinearitas
#Multicollinearity
vif(model_ols)
## X1 X2 X3
## 1.406373 1.184105 1.210866
Uji Autokorelasi
CoordK <- coordinates(Jatim)
#KNN=5
W.K_5 <- knn2nb(knearneigh(CoordK, k = 5 ), row.names = Id)
WL.K_5 <- nb2listw(W.K_5, style = 'W') ; WL.K_5
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 190
## Percentage nonzero weights: 13.15789
## Average number of links: 5
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 13.44 156.96
#KNN=6
W.K_6 <- knn2nb(knearneigh(CoordK, k = 6 ), row.names = Id)
WL.K_6 <- nb2listw(W.K_6, style = 'W') ; WL.K_6
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 228
## Percentage nonzero weights: 15.78947
## Average number of links: 6
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 38 1444 38 11.22222 156.8333
#Rook
W.R <- poly2nb(Jatim, row.names = Id, queen = FALSE)
## Warning in poly2nb(Jatim, row.names = Id, queen = FALSE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(Jatim, row.names = Id, queen = FALSE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
WL.R <- nb2listw(W.R, style = 'W', zero.policy = TRUE) ; WL.R
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 78
## Percentage nonzero weights: 5.401662
## Average number of links: 2.052632
## 6 regions with no links:
## 11, 12, 14, 15, 17, 20
## 10 disjoint connected subgraphs
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 32 1024 32 29.83333 141.3244
#Queen
W.Q <- poly2nb(Jatim, row.names = Id, queen = TRUE)
## Warning in poly2nb(Jatim, row.names = Id, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in poly2nb(Jatim, row.names = Id, queen = TRUE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
WL.Q <- nb2listw(W.Q, style = 'W', zero.policy = TRUE) ; WL.Q
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 38
## Number of nonzero links: 78
## Percentage nonzero weights: 5.401662
## Average number of links: 2.052632
## 6 regions with no links:
## 11, 12, 14, 15, 17, 20
## 10 disjoint connected subgraphs
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 32 1024 32 29.83333 141.3244
#Choosing spatial weights
MI.knn5 <- moran.test(residuals, WL.K_5)
MI.knn6 <- moran.test(residuals, WL.K_6)
MI.rook <- moran.test(residuals, WL.R)
MI.queen <- moran.test(residuals, WL.Q)
MI.knn5$estimate
## Moran I statistic Expectation Variance
## 0.040160884 -0.027027027 0.007936345
MI.knn6$estimate
## Moran I statistic Expectation Variance
## 0.023152184 -0.027027027 0.006374622
MI.rook$estimate
## Moran I statistic Expectation Variance
## 0.23768940 -0.03225806 0.02727792
MI.queen$estimate
## Moran I statistic Expectation Variance
## 0.23768940 -0.03225806 0.02727792
moran_index<-data.frame("Weight matrix"=c("KNN (k = 5)", "KNN (k = 6)", "Rook Contiguity", "Queen Contiguity"), "Moran's I"=c(0.040160884, 0.023152184, 0.23768940, 0.23768940))
moran_index
## Weight.matrix Moran.s.I
## 1 KNN (k = 5) 0.04016088
## 2 KNN (k = 6) 0.02315218
## 3 Rook Contiguity 0.23768940
## 4 Queen Contiguity 0.23768940
W_optimum <- WL.Q
moran.test(Y,W_optimum)
##
## Moran I test under randomisation
##
## data: Y
## weights: W_optimum
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 1.9973, p-value = 0.0229
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.26798372 -0.03225806 0.02259684
Uji Heteroskedastisitas
# Heteroskedasticity test (Breusch-Pagan)
bptest(model_ols)
##
## studentized Breusch-Pagan test
##
## data: model_ols
## BP = 7.8646, df = 3, p-value = 0.04889
Model GWR
Kernel Adaptive Gaussian
#Kernel Adaptive Gaussian
adapt_bandwidth <- gwr.sel(Y ~ X1 + X2 + X3, data = Jatim, coords = coords, adapt=TRUE)
## Warning in gwr.sel(Y ~ X1 + X2 + X3, data = Jatim, coords = coords, adapt =
## TRUE): data is Spatial* object, ignoring coords argument
## Adaptive q: 0.381966 CV score: 3039334
## Adaptive q: 0.618034 CV score: 3161403
## Adaptive q: 0.236068 CV score: 2975060
## Adaptive q: 0.145898 CV score: 2969774
## Adaptive q: 0.1728634 CV score: 2951236
## Adaptive q: 0.1884997 CV score: 2953267
## Adaptive q: 0.1772967 CV score: 2951639
## Adaptive q: 0.1625635 CV score: 2955951
## Adaptive q: 0.1689292 CV score: 2951850
## Adaptive q: 0.173541 CV score: 2951233
## Adaptive q: 0.1732817 CV score: 2951231
## Adaptive q: 0.173241 CV score: 2951231
## Adaptive q: 0.1733224 CV score: 2951231
## Adaptive q: 0.1732817 CV score: 2951231
adapt_bandwidth
## [1] 0.1732817
gwr.model1 <- gwr(Y ~ X1 + X2 + X3,
data=Jatim,
adapt=adapt_bandwidth,
hatmatrix = TRUE,
se.fit = TRUE)
print(gwr.model1)
## Call:
## gwr(formula = Y ~ X1 + X2 + X3, data = Jatim, adapt = adapt_bandwidth,
## hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Adaptive quantile: 0.1732817 (about 6 of 38 data points)
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. -1317.3316 -839.8762 -622.3795 -458.3273 -240.6965 -551.0948
## X1 -11.2238 -8.8537 -6.6590 -5.1017 -1.8238 -8.7067
## X2 33.8639 40.8100 48.3670 51.9116 54.5826 48.5485
## X3 133.7041 159.2907 198.2168 207.2017 214.3086 179.0434
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 14.79861
## Effective degrees of freedom (residual: 2traceS - traceS'S): 23.20139
## Sigma (residual: 2traceS - traceS'S): 244.2271
## Effective number of parameters (model: traceS): 11.19937
## Effective degrees of freedom (model: traceS): 26.80063
## Sigma (model: traceS): 227.2366
## Sigma (ML): 190.8354
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 544.3308
## AIC (GWR p. 96, eq. 4.22): 518.146
## Residual sum of squares: 1383890
## Quasi-global R2: 0.9395351
Kernel Fixed Gaussian
#Kernel Fixed Gaussian
bwl <- gwr.sel(Y ~ X1 + X2 + X3, data=Jatim, coords = coords)
## Warning in gwr.sel(Y ~ X1 + X2 + X3, data = Jatim, coords = coords): data is
## Spatial* object, ignoring coords argument
## Bandwidth: 141.0182 CV score: 3077925
## Bandwidth: 227.9445 CV score: 3293379
## Bandwidth: 87.29485 CV score: 2944328
## Bandwidth: 54.09198 CV score: 3210295
## Bandwidth: 103.8601 CV score: 2970894
## Bandwidth: 74.61248 CV score: 2943064
## Bandwidth: 79.98466 CV score: 2940415
## Bandwidth: 80.33917 CV score: 2940443
## Bandwidth: 79.76115 CV score: 2940407
## Bandwidth: 79.66706 CV score: 2940406
## Bandwidth: 79.67186 CV score: 2940406
## Bandwidth: 79.67232 CV score: 2940406
## Bandwidth: 79.67228 CV score: 2940406
## Bandwidth: 79.67237 CV score: 2940406
## Bandwidth: 79.67232 CV score: 2940406
bwl
## [1] 79.67232
gwr.model2 <- gwr(Y ~ X1 + X2 + X3,
data=Jatim,
bandwidth=bwl,
hatmatrix = TRUE,
se.fit = TRUE)
print(gwr.model2)
## Call:
## gwr(formula = Y ~ X1 + X2 + X3, data = Jatim, bandwidth = bwl,
## hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 79.67232
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. -1021.7500 -759.8757 -662.4562 -597.5695 -13.3988 -551.0948
## X1 -12.9659 -8.4205 -7.3799 -6.7843 -3.1205 -8.7067
## X2 37.9597 44.6607 48.0841 50.5843 56.2841 48.5485
## X3 111.8889 162.6385 193.3735 210.4413 214.2912 179.0434
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 11.47622
## Effective degrees of freedom (residual: 2traceS - traceS'S): 26.52378
## Sigma (residual: 2traceS - traceS'S): 248.2305
## Effective number of parameters (model: traceS): 8.868861
## Effective degrees of freedom (model: traceS): 29.13114
## Sigma (model: traceS): 236.8614
## Sigma (ML): 207.3868
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 540.9126
## AIC (GWR p. 96, eq. 4.22): 522.1367
## Residual sum of squares: 1634353
## Quasi-global R2: 0.9285918
Kernel Fixed Bi-square
#Bi-square
bw2 <- gwr.sel(Y ~ X1 + X2 + X3, data = Jatim, coords = coords, gweight = gwr.bisquare)
## Warning in gwr.sel(Y ~ X1 + X2 + X3, data = Jatim, coords = coords, gweight =
## gwr.bisquare): data is Spatial* object, ignoring coords argument
## Bandwidth: 141.0182 CV score: 3610540
## Bandwidth: 227.9445 CV score: 2931008
## Bandwidth: 281.6679 CV score: 2978615
## Bandwidth: 247.6461 CV score: 2944700
## Bandwidth: 194.7416 CV score: 2897358
## Bandwidth: 174.2211 CV score: 2971086
## Bandwidth: 207.424 CV score: 2906377
## Bandwidth: 186.9035 CV score: 2906735
## Bandwidth: 197.2574 CV score: 2897711
## Bandwidth: 194.7508 CV score: 2897358
## Bandwidth: 195.1714 CV score: 2897344
## Bandwidth: 195.9682 CV score: 2897403
## Bandwidth: 195.1413 CV score: 2897344
## Bandwidth: 195.1437 CV score: 2897344
## Bandwidth: 195.1436 CV score: 2897344
## Bandwidth: 195.1435 CV score: 2897344
## Bandwidth: 195.1436 CV score: 2897344
## Bandwidth: 195.1436 CV score: 2897344
bw2
## [1] 195.1436
gwr.model3 <- gwr(Y ~ X1 + X2 + X3,
data = Jatim,
bandwidth=bw2,
gweight = gwr.bisquare,
hatmatrix = TRUE,
se.fit = TRUE)
print(gwr.model3)
## Call:
## gwr(formula = Y ~ X1 + X2 + X3, data = Jatim, bandwidth = bw2,
## gweight = gwr.bisquare, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.bisquare
## Fixed bandwidth: 195.1436
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. -1080.2463 -799.5010 -691.4953 -596.5033 515.9442 -551.0948
## X1 -17.3820 -8.1156 -6.9983 -6.3985 -3.4894 -8.7067
## X2 33.2715 45.1390 47.9548 50.6643 57.5917 48.5485
## X3 96.2337 159.0229 194.4323 212.0990 217.8687 179.0434
## Number of data points: 38
## Effective number of parameters (residual: 2traceS - traceS'S): 10.73944
## Effective degrees of freedom (residual: 2traceS - traceS'S): 27.26056
## Sigma (residual: 2traceS - traceS'S): 246.9093
## Effective number of parameters (model: traceS): 8.639503
## Effective degrees of freedom (model: traceS): 29.3605
## Sigma (model: traceS): 237.9157
## Sigma (ML): 209.1284
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 540.6793
## AIC (GWR p. 96, eq. 4.22): 522.5429
## Residual sum of squares: 1661918
## Quasi-global R2: 0.9273875
AIC <- c(gwr.model1$results$AICh, gwr.model2$results$AICh, gwr.model3$results$AICh)
AIC
## [1] 518.1460 522.1367 522.5429
#Model yg dipih (model 1)
hasil <- as.data.frame(gwr.model1$SDF)
head(hasil)
## sum.w X.Intercept. X1 X2 X3 X.Intercept._se X1_se
## 1 12.265198 -730.1325 -5.407390 53.85055 153.2463 279.9975 3.971981
## 2 15.004220 -431.3569 -8.432239 50.73394 156.5619 260.3502 3.553628
## 3 11.244170 -819.7430 -6.664770 47.29370 205.4431 396.6287 5.134289
## 4 9.274121 -449.9407 -10.516800 40.68505 206.3715 471.4948 5.661476
## 5 11.155135 -1073.5792 -3.907525 45.49460 209.2322 411.3830 5.155670
## 6 12.316979 -446.5919 -7.315188 51.96379 144.8202 274.8444 3.951762
## X2_se X3_se gwr.e pred pred.se localR2
## 1 3.857264 24.78224 -481.61241 1406.6124 83.07099 0.9346639
## 2 3.534635 22.25988 -87.38474 2086.3847 83.26131 0.9374351
## 3 3.996867 27.54708 34.13294 122.8671 88.48668 0.9505882
## 4 5.273913 31.40162 101.84083 521.1592 104.08301 0.9462961
## 5 4.329789 25.64478 61.87939 1368.1206 75.69445 0.9287710
## 6 3.814412 24.46196 -74.03563 927.0356 71.19907 0.9337385
## X.Intercept._se_EDF X1_se_EDF X2_se_EDF X3_se_EDF pred.se.1
## 1 300.9329 4.268966 4.145672 26.63520 89.28221
## 2 279.8166 3.819332 3.798920 23.92426 89.48676
## 3 426.2846 5.518179 4.295713 29.60677 95.10283
## 4 506.7485 6.084785 5.668243 33.74952 111.86530
## 5 442.1421 5.541159 4.653527 27.56224 81.35413
## 6 295.3945 4.247235 4.099616 26.29098 76.52262
# Load packages
library(spgwr)
library(tmap)
## Warning: package 'tmap' was built under R version 4.3.3
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(sp)
# Gabungkan hasil GWR dengan data spatial
gwr.map <- cbind(Jatim, as.matrix(hasil))
tm_shape(gwr.map) +
tm_polygons("localR2", title = "Local R²")
# Plot residuals dari GWR model 1
resid_gwr <- hasil$gwr.e # Pastikan `gwr.e` ada di `hasil`
Jatim$resid_gwr <- resid_gwr
spplot(Jatim, zcol = "resid_gwr", main = "Peta sebaran residual model GWR")
Uji Normalitas
#Asumsi
#Normalitas pada model GWR
shapiro.test(resid_gwr)
##
## Shapiro-Wilk normality test
##
## data: resid_gwr
## W = 0.92587, p-value = 0.0149
# Uji Heteroskedastisitas Breusch-Pagan pada model GWR
bp_test <- bptest(gwr.model1$SDF$gwr.e ~ gwr.model1$SDF$pred +
I(gwr.model1$SDF$pred^2))
# Menampilkan hasil uji
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: gwr.model1$SDF$gwr.e ~ gwr.model1$SDF$pred + I(gwr.model1$SDF$pred^2)
## BP = 8.2264, df = 2, p-value = 0.01636
#Choosing spatial weights
MI.knn5 <- moran.test(resid_gwr, WL.K_5)
MI.knn6 <- moran.test(resid_gwr, WL.K_6)
MI.rook1 <- moran.test(resid_gwr, WL.R)
MI.queen1 <- moran.test(resid_gwr, WL.Q)
MI.knn5$estimate
## Moran I statistic Expectation Variance
## -0.140208956 -0.027027027 0.007664306
MI.knn6$estimate
## Moran I statistic Expectation Variance
## -0.150617495 -0.027027027 0.006156682
MI.rook1$estimate
## Moran I statistic Expectation Variance
## -0.04494032 -0.03225806 0.02614554
MI.queen1$estimate
## Moran I statistic Expectation Variance
## -0.04494032 -0.03225806 0.02614554
W_optimum <- WL.Q
gwr.morantest(gwr.model1, W_optimum, zero.policy = TRUE)
##
## Leung et al. 2000 three moment approximation for Moran's I
##
## data: GWR residuals
## statistic = 679.33, df = 694.06, p-value = 0.3519
## sample estimates:
## I
## -0.04491921
moran.test(resid_gwr,W_optimum)
##
## Moran I test under randomisation
##
## data: resid_gwr
## weights: W_optimum
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = -0.078433, p-value = 0.5313
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04494032 -0.03225806 0.02614554
#Comparing matrix
data.frame("Model" = c("GWR", "OLS"),
"AIC" = c(gwr.model1$results$AICh, AIC(model_ols)),
"R2" = c(0.9395351, 0.8742)) %>% arrange(AIC)
## Model AIC R2
## 1 GWR 518.1460 0.9395351
## 2 OLS 541.5679 0.8742000
##Membuat plot -gglot2##
gwr_df <- as.data.frame(gwr.model1$SDF)
gwr_df
## sum.w X.Intercept. X1 X2 X3 X.Intercept._se X1_se
## 1 12.265198 -730.1325 -5.407390 53.85055 153.2463 279.9975 3.971981
## 2 15.004220 -431.3569 -8.432239 50.73394 156.5619 260.3502 3.553628
## 3 11.244170 -819.7430 -6.664770 47.29370 205.4431 396.6287 5.134289
## 4 9.274121 -449.9407 -10.516800 40.68505 206.3715 471.4948 5.661476
## 5 11.155135 -1073.5792 -3.907525 45.49460 209.2322 411.3830 5.155670
## 6 12.316979 -446.5919 -7.315188 51.96379 144.8202 274.8444 3.951762
## 7 9.668275 -1022.3635 -4.608897 52.71042 186.0124 290.3973 4.285834
## 8 11.952071 -415.8806 -8.369501 51.31549 154.3498 277.1473 3.962454
## 9 11.500410 -1175.9066 -3.036596 46.81401 214.1104 396.3312 5.058747
## 10 11.116893 -636.4055 -8.380267 40.16436 210.5396 462.8341 5.490538
## 11 11.564701 -526.4166 -10.054712 41.18484 211.2899 432.3451 5.277304
## 12 10.855918 -552.2195 -8.954318 38.76984 207.3243 482.1185 5.665782
## 13 10.750124 -475.4134 -9.493278 37.22645 205.6913 509.7459 5.976471
## 14 11.918693 -676.4753 -7.585958 47.72301 194.4045 372.7648 4.871255
## 15 10.013050 -1317.3316 -1.823796 49.59618 213.8035 392.1906 5.426254
## 16 10.612038 -754.0418 -6.329215 51.75515 182.2389 290.6864 4.306759
## 17 13.977832 -611.9353 -6.354265 52.26082 159.1226 262.9225 3.702177
## 18 10.448523 -1217.4614 -2.155869 50.94154 193.7140 326.4834 4.455480
## 19 9.327570 -564.7655 -6.619877 51.68720 159.7948 325.1014 4.503312
## 20 9.398477 -332.0094 -10.312369 34.99514 198.9166 562.9513 6.339120
## 21 10.663300 -492.2885 -9.472982 37.68391 206.8338 503.5698 5.915442
## 22 11.898324 -582.4619 -8.551877 47.11272 192.2857 367.6267 4.769616
## 23 11.542353 -1109.5288 -4.013111 49.01104 210.6925 355.9214 4.921974
## 24 13.039184 -846.5873 -6.467794 42.39529 214.3086 416.0818 5.066514
## 25 11.713266 -744.0512 -7.466721 41.49246 211.7221 426.5000 5.284731
## 26 9.761980 -424.6708 -10.402799 37.58986 207.4579 510.2258 5.770209
## 27 13.318259 -632.8237 -4.660788 54.12848 134.2273 273.5117 3.964517
## 28 10.376510 -774.1076 -6.210819 50.93491 187.4365 331.6820 4.755038
## 29 9.572633 -240.6965 -11.223779 33.86393 197.5169 575.3969 6.427969
## 30 13.618246 -565.9372 -6.653304 52.25051 155.8174 262.7647 3.714947
## 31 12.478250 -673.9679 -4.600556 54.58260 137.1437 280.0302 4.051704
## 32 9.965157 -913.8225 -6.261327 51.38720 200.2493 300.0328 4.596569
## 33 11.114704 -452.6320 -6.814412 52.34322 139.0854 280.3911 4.086306
## 34 13.312548 -592.2525 -4.999780 53.54905 133.7041 268.8090 3.874869
## 35 10.118450 -849.0901 -6.529179 52.85175 184.9123 283.2260 4.249143
## 36 10.302410 -339.8493 -11.099648 36.47770 204.7362 518.5812 5.916990
## 37 13.483308 -1016.8865 -4.616926 47.71115 201.3304 306.9940 3.991820
## 38 10.678984 -383.9151 -10.916572 37.62634 206.5942 486.3062 5.726127
## X2_se X3_se gwr.e pred pred.se localR2
## 1 3.857264 24.78224 -481.612408 1406.6124 83.07099 0.9346639
## 2 3.534635 22.25988 -87.384742 2086.3847 83.26131 0.9374351
## 3 3.996867 27.54708 34.132935 122.8671 88.48668 0.9505882
## 4 5.273913 31.40162 101.840827 521.1592 104.08301 0.9462961
## 5 4.329789 25.64478 61.879394 1368.1206 75.69445 0.9287710
## 6 3.814412 24.46196 -74.035626 927.0356 71.19907 0.9337385
## 7 4.315700 30.14368 -406.673231 1843.6732 71.54707 0.9378409
## 8 3.795377 24.80459 208.615357 2553.3846 114.15726 0.9365325
## 9 4.034937 26.26222 -440.273665 1726.2737 72.26431 0.9381461
## 10 4.914013 26.84870 -47.599823 1466.5998 109.76690 0.9294031
## 11 4.823009 27.18907 -67.954468 301.9545 106.24078 0.9435889
## 12 5.199493 27.41296 128.284445 440.7156 110.92551 0.9231082
## 13 5.651059 33.75132 34.385770 491.6142 146.42732 0.9195462
## 14 3.865646 26.77598 -3.573223 1363.5732 120.01459 0.9469732
## 15 4.235637 29.64770 -7.142848 329.1428 97.21675 0.9449966
## 16 4.029384 28.04841 175.612397 310.3876 83.20564 0.9468749
## 17 3.622499 23.83647 -104.689525 370.6895 89.82804 0.9362270
## 18 4.133026 28.54932 222.510737 1270.4893 90.71495 0.9339739
## 19 4.027355 28.88973 315.013586 802.9864 98.38942 0.9327812
## 20 6.267941 38.35175 -178.082270 768.0823 80.18112 0.9021339
## 21 5.580622 33.69091 43.272610 415.7274 91.22356 0.9271630
## 22 3.834018 26.56139 -33.079818 1840.0798 119.45107 0.9448456
## 23 4.004173 28.11722 -85.823077 1140.8231 96.16301 0.9475888
## 24 4.415833 24.75643 129.222845 587.7772 57.68380 0.9291290
## 25 4.672389 28.07602 -39.174616 795.1746 54.55730 0.9343443
## 26 5.592429 33.34668 39.571923 247.4281 131.75637 0.9366475
## 27 3.560254 24.47204 123.055700 598.9443 84.05302 0.9285142
## 28 4.076983 28.63460 182.026597 1577.9734 76.94622 0.9486137
## 29 6.609263 39.16794 -26.726258 925.7263 88.63606 0.9158415
## 30 3.616443 23.91725 -386.218571 1538.2186 56.58898 0.9350866
## 31 3.706137 24.70814 -35.869852 850.8699 181.94820 0.9305039
## 32 4.231130 30.94431 257.873283 2263.1267 151.41392 0.9508592
## 33 3.914967 25.15022 176.589818 789.4102 77.77259 0.9316765
## 34 3.476168 24.26889 256.653958 1080.3460 88.99932 0.9277630
## 35 4.254075 29.57888 241.407903 3859.5921 172.33537 0.9440987
## 36 5.912719 32.04724 -25.932092 426.9321 82.51082 0.9314526
## 37 3.746241 23.20926 -68.153488 1324.1535 80.97229 0.9358330
## 38 5.647129 28.55261 -130.939976 957.9400 113.57017 0.9332491
## X.Intercept._se_EDF X1_se_EDF X2_se_EDF X3_se_EDF pred.se.1
## 1 300.9329 4.268966 4.145672 26.63520 89.28221
## 2 279.8166 3.819332 3.798920 23.92426 89.48676
## 3 426.2846 5.518179 4.295713 29.60677 95.10283
## 4 506.7485 6.084785 5.668243 33.74952 111.86530
## 5 442.1421 5.541159 4.653527 27.56224 81.35413
## 6 295.3945 4.247235 4.099616 26.29098 76.52262
## 7 312.1103 4.606286 4.638385 32.39753 76.89665
## 8 297.8696 4.258727 4.079157 26.65923 122.69280
## 9 425.9649 5.436989 4.336629 28.22585 77.66752
## 10 497.4402 5.901066 5.281434 28.85618 117.97417
## 11 464.6715 5.671888 5.183625 29.22200 114.18441
## 12 518.1665 6.089412 5.588259 29.46263 119.21941
## 13 547.8596 6.423332 6.073588 36.27491 157.37569
## 14 400.6365 5.235479 4.154680 28.77801 128.98808
## 15 421.5147 5.831975 4.552335 31.86445 104.48565
## 16 312.4210 4.628775 4.330661 30.14559 89.42692
## 17 282.5812 3.978988 3.893354 25.61872 96.54449
## 18 350.8946 4.788616 4.442052 30.68395 97.49771
## 19 349.4093 4.840025 4.328480 31.04981 105.74600
## 20 605.0432 6.813096 6.736594 41.21931 86.17626
## 21 541.2217 6.357739 5.997885 36.20998 98.04434
## 22 395.1141 5.126241 4.120688 28.54739 128.38243
## 23 382.5336 5.289990 4.303565 30.21954 103.35312
## 24 447.1922 5.445337 4.746005 26.60747 61.99682
## 25 458.3894 5.679871 5.021744 30.17526 58.63655
## 26 548.3754 6.201648 6.010574 35.84001 141.60779
## 27 293.9621 4.260944 3.826454 26.30181 90.33766
## 28 356.4819 5.110572 4.381819 30.77561 82.69949
## 29 618.4193 6.908588 7.103438 42.09653 95.26338
## 30 282.4116 3.992714 3.886845 25.70555 60.82014
## 31 300.9681 4.354650 3.983245 26.55557 195.55246
## 32 322.4662 4.940255 4.547492 33.25802 162.73513
## 33 301.3559 4.391839 4.207690 27.03070 83.58765
## 34 288.9078 4.164592 3.736081 26.08347 95.65379
## 35 304.4028 4.566851 4.572152 31.79049 185.22088
## 36 557.3555 6.359404 6.354813 34.44341 88.68016
## 37 329.9480 4.290288 4.026347 24.94462 87.02659
## 38 522.6673 6.154270 6.069365 30.68749 122.06181
Jatim_merged1 <- Jatim_merged %>%
bind_cols(gwr_df %>%
select(X.Intercept., X1, X2, X3))
# Membuat peta dengan geom_sf
ggplot(data = Jatim_merged1) +
geom_sf(aes(fill = X.Intercept.), color = "blue") +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabptn..Kt), size = 2, color = "white") +
labs(title = "Koefisien GWR untuk intercept", fill = "Koefisien Beta0") +
theme_minimal()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# Membuat peta dengan geom_sf
ggplot(data = Jatim_merged1) +
geom_sf(aes(fill = X1), color = "blue") +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabptn..Kt), size = 2, color = "white") +
labs(title = "Koefisien GWR untuk X1", fill = "Koefisien Beta1") +
theme_minimal()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# Membuat peta dengan geom_sf
ggplot(data = Jatim_merged1) +
geom_sf(aes(fill = X2), color = "blue") +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabptn..Kt), size = 2, color = "white") +
labs(title = "Koefisien GWR untuk X2", fill = "Koefisien Beta2") +
theme_minimal()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
# Membuat peta dengan geom_sf
ggplot(data = Jatim_merged1) +
geom_sf(aes(fill = X3), color = "blue") +
scale_fill_viridis_c() +
geom_sf_text(aes(label = Kabptn..Kt), size = 2, color = "white") +
labs(title = "Koefisien GWR untuk X3", fill = "Koefisien Beta3") +
theme_minimal()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
#Local significance test
t1 <- gwr.model1$SDF$X1/gwr.model1$SDF$X1_se
t2 <- gwr.model2$SDF$X2/gwr.model2$SDF$X2_se
t3 <- gwr.model3$SDF$X3/gwr.model3$SDF$X3_se
pvalue1 <- 2*pt(q=t1, df=38-1, lower.tail = FALSE)
pvalue2 <- 2*pt(q=t2, df=38-1, lower.tail = FALSE)
pvalue3 <- 2*pt(q=t3, df=38-1, lower.tail = FALSE)
#Variabel X1: Persentase TTU yang memenuhi syarat kesehatan
# Assuming t1 and pvalue1 are already defined and have 38 elements each
Id <- seq_len(length(t1)) # Create Id with 38 elements
# Create data frame
data.frame(Id, t1, pvalue1) %>%
mutate(ket = ifelse(pvalue1 < 0.05, "Reject", "Not reject"))
## Id t1 pvalue1 ket
## 1 1 -1.3613837 1.818374 Not reject
## 2 2 -2.3728540 1.977040 Not reject
## 3 3 -1.2980903 1.797710 Not reject
## 4 4 -1.8576073 1.928803 Not reject
## 5 5 -0.7579083 1.546694 Not reject
## 6 6 -1.8511205 1.927853 Not reject
## 7 7 -1.0753792 1.710832 Not reject
## 8 8 -2.1122015 1.958519 Not reject
## 9 9 -0.6002664 1.448012 Not reject
## 10 10 -1.5263107 1.864564 Not reject
## 11 11 -1.9052744 1.935462 Not reject
## 12 12 -1.5804206 1.877477 Not reject
## 13 13 -1.5884420 1.879304 Not reject
## 14 14 -1.5572904 1.872085 Not reject
## 15 15 -0.3361060 1.261310 Not reject
## 16 16 -1.4696007 1.849878 Not reject
## 17 17 -1.7163592 1.905541 Not reject
## 18 18 -0.4838692 1.368670 Not reject
## 19 19 -1.4700019 1.849986 Not reject
## 20 20 -1.6267824 1.887729 Not reject
## 21 21 -1.6013988 1.882207 Not reject
## 22 22 -1.7929905 1.918847 Not reject
## 23 23 -0.8153459 1.579908 Not reject
## 24 24 -1.2765766 1.790294 Not reject
## 25 25 -1.4128857 1.833953 Not reject
## 26 26 -1.8028460 1.920438 Not reject
## 27 27 -1.1756256 1.752746 Not reject
## 28 28 -1.3061555 1.800438 Not reject
## 29 29 -1.7460848 1.910905 Not reject
## 30 30 -1.7909553 1.918516 Not reject
## 31 31 -1.1354620 1.736516 Not reject
## 32 32 -1.3621739 1.818621 Not reject
## 33 33 -1.6676214 1.896166 Not reject
## 34 34 -1.2903095 1.795051 Not reject
## 35 35 -1.5365874 1.867097 Not reject
## 36 36 -1.8758943 1.931424 Not reject
## 37 37 -1.1565969 1.745150 Not reject
## 38 38 -1.9064495 1.935619 Not reject
#Variabel X2: Jumlah Puskesmas
# Assuming t2 and pvalue2 are already defined and have 38 elements each
Id <- seq_len(length(t2)) # Create Id with 38 elements
# Create data frame
data.frame(Id, t2, pvalue2) %>%
mutate(ket = ifelse(pvalue2 < 0.05, "Reject", "Not reject"))
## Id t2 pvalue2 ket
## 1 1 13.873848 3.081711e-16 Reject
## 2 2 8.017502 1.309017e-09 Reject
## 3 3 13.351316 1.012920e-15 Reject
## 4 4 11.609183 6.755123e-14 Reject
## 5 5 11.766824 4.549088e-14 Reject
## 6 6 11.226486 1.787200e-13 Reject
## 7 7 13.459596 7.895506e-16 Reject
## 8 8 11.974800 2.713078e-14 Reject
## 9 9 12.914548 2.805427e-15 Reject
## 10 10 11.790911 4.283590e-14 Reject
## 11 11 11.395665 1.159796e-13 Reject
## 12 12 11.488706 9.157581e-14 Reject
## 13 13 9.747504 9.169703e-12 Reject
## 14 14 13.469495 7.718233e-16 Reject
## 15 15 13.371581 9.666790e-16 Reject
## 16 16 14.169532 1.593087e-16 Reject
## 17 17 14.327249 1.124845e-16 Reject
## 18 18 13.092090 1.849225e-15 Reject
## 19 19 14.072189 1.977491e-16 Reject
## 20 20 10.339597 1.831759e-12 Reject
## 21 21 8.768660 1.451421e-10 Reject
## 22 22 13.260117 1.250715e-15 Reject
## 23 23 13.447593 8.116054e-16 Reject
## 24 24 11.717295 5.149052e-14 Reject
## 25 25 9.234997 3.834963e-11 Reject
## 26 26 6.667853 7.884708e-08 Reject
## 27 27 14.213821 1.444366e-16 Reject
## 28 28 14.055490 2.052412e-16 Reject
## 29 29 8.731104 1.617542e-10 Reject
## 30 30 14.167671 1.599669e-16 Reject
## 31 31 14.172820 1.581525e-16 Reject
## 32 32 13.838853 3.334143e-16 Reject
## 33 33 11.097279 2.492688e-13 Reject
## 34 34 13.736058 4.204996e-16 Reject
## 35 35 13.794754 3.682622e-16 Reject
## 36 36 8.630168 2.166268e-10 Reject
## 37 37 12.197957 1.567520e-14 Reject
## 38 38 10.023718 4.301027e-12 Reject
#Variabel x3: Tingkat Pengangguran Terbuka
# Assuming t3 and pvalue3 are already defined and have 38 elements each
Id <- seq_len(length(t3)) # Create Id with 38 elements
# Create data frame
data.frame(Id, t3, pvalue3) %>%
mutate(ket = ifelse(pvalue3 < 0.05, "Reject", "Not reject"))
## Id t3 pvalue3 ket
## 1 1 6.058561 5.231051e-07 Reject
## 2 2 1.582605 1.220231e-01 Not reject
## 3 3 8.524360 2.946225e-10 Reject
## 4 4 9.010140 7.261022e-11 Reject
## 5 5 9.014950 7.162063e-11 Reject
## 6 6 4.499305 6.545054e-05 Reject
## 7 7 7.617543 4.332451e-09 Reject
## 8 8 5.172453 8.269736e-06 Reject
## 9 9 8.946571 8.707156e-11 Reject
## 10 10 9.267264 3.501112e-11 Reject
## 11 11 9.098794 5.641072e-11 Reject
## 12 12 9.291240 3.272299e-11 Reject
## 13 13 8.470866 3.443568e-10 Reject
## 14 14 8.246495 6.648757e-10 Reject
## 15 15 8.490406 3.252722e-10 Reject
## 16 16 7.245551 1.337787e-08 Reject
## 17 17 6.461136 1.495543e-07 Reject
## 18 18 8.230508 6.969426e-10 Reject
## 19 19 6.730281 6.501857e-08 Reject
## 20 20 8.916864 9.480078e-11 Reject
## 21 21 7.306582 1.110893e-08 Reject
## 22 22 8.177088 8.159079e-10 Reject
## 23 23 8.455625 3.600287e-10 Reject
## 24 24 9.291283 3.271895e-11 Reject
## 25 25 7.391838 8.573576e-09 Reject
## 26 26 4.291311 1.226221e-04 Reject
## 27 27 4.580703 5.110192e-05 Reject
## 28 28 7.571718 4.974442e-09 Reject
## 29 29 7.811524 2.419490e-09 Reject
## 30 30 6.213707 3.226501e-07 Reject
## 31 31 5.135561 9.272069e-06 Reject
## 32 32 7.680269 3.587051e-09 Reject
## 33 33 4.078105 2.315378e-04 Reject
## 34 34 3.652535 7.990541e-04 Reject
## 35 35 7.308961 1.102881e-08 Reject
## 36 36 7.900472 1.854710e-09 Reject
## 37 37 8.331196 5.182934e-10 Reject
## 38 38 8.949932 8.623836e-11 Reject