INPUT DATA DAN PLOT PETA

setwd("D:/KULIAH/Spasial")
st_layers("gadm41_USA.gpkg")

## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields crs_name
## 1  ADM_ADM_0 Multi Polygon        1      2   WGS 84
## 2  ADM_ADM_1 Multi Polygon       51     11   WGS 84
## 3  ADM_ADM_2 Multi Polygon     3148     13   WGS 84

USA = st_read("gadm41_USA.gpkg", layer = "ADM_ADM_2")

## Reading layer `ADM_ADM_2' from data source `D:\KULIAH\Spasial\gadm41_USA.gpkg' using driver `GPKG'
## Simple feature collection with 3148 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1506 ymin: 18.90986 xmax: 179.7734 ymax: 72.6875
## Geodetic CRS:  WGS 84

Miss <- USA[USA$NAME_1 == "Mississippi", ]
centroid <- st_centroid(Miss)

## Warning: st_centroid assumes attributes are constant over geometries

Coordk <- st_coordinates(centroid)
plot(st_geometry(Miss), axes=TRUE, col="gray90")
text(Coordk[,1],Coordk[,2],Miss$NAME_2, col="black",cex=0.42)

data <- read.csv("Mississippi.csv")
head(data)

##    County Poverty.Rate..y. Unemployed.Rate....
## 1   Adams             27.2                 7.1
## 2  Alcorn             17.1                 6.4
## 3   Amite             27.1                 7.5
## 4  Attala             23.6                 7.2
## 5  Benton             17.6                 4.4
## 6 Bolivar             31.8                 7.4
##   Median.household.income..dollars. Education.less.than.highschool....
## 1                             37271                               17.5
## 2                             47716                               18.1
## 3                             34866                               12.5
## 4                             42680                               15.9
## 5                             38750                               22.0
## 6                             37845                               17.8

PLOT PERSEBARAN TIAP VARIABEL

Miss$poverty = data$Poverty.Rate..y.
Miss$income = data$Median.household.income..dollars.
Miss$unemployment = data$Unemployed.Rate....
Miss$education = data$Education.less.than.highschool....

1. Persebaran Kemiskinan

tm_shape(Miss) + tm_polygons("poverty",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("NAME_2",size=0.5)

2. Persebaran Penghasilan Rumah Tangga

tm_shape(Miss) + tm_polygons("income",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("NAME_2",size=0.5)

3. Persebaran Pengangguran

tm_shape(Miss) + tm_polygons("unemployment",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("NAME_2",size=0.5)

4. Persebaran Edukasi dibawah SMA

tm_shape(Miss) + tm_polygons("education",style="cont")+ tm_borders(lwd=1,col="black") + tm_text("NAME_2",size=0.5)

MEMBANGUN MATRIKS PEMBOBOT SPASIAL DENGAN ROOK

Mendapatkan Koordinat Centroid

centroid <- st_centroid(Miss)

## Warning: st_centroid assumes attributes are constant over geometries

Coordk <- st_coordinates(centroid)

Mendapatkan Matriks Pembobot Spasial

W = poly2nb(Miss, row.names = Miss$ID, queen = FALSE)
WB = nb2mat(W, style = 'B', zero.policy = TRUE)
WB[1:5,1:5]

##      [,1] [,2] [,3] [,4] [,5]
## 1405    0    0    0    0    0
## 1406    0    0    0    0    0
## 1407    0    0    0    0    0
## 1408    0    0    0    0    0
## 1409    0    0    0    0    0

WL = nb2listw(W)
WL

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 82 
## Number of nonzero links: 404 
## Percentage nonzero weights: 6.008328 
## Average number of links: 4.926829 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 82 6724 82 34.95435 332.2967

Plot Hubungan Lokasi pada Peta

plot(st_geometry(Miss), axes=TRUE, col="gray90")
plot(W,Coordk,col="red",add=TRUE)

Moran’s Test

poverty = data$Poverty.Rate..y.
moran.test(poverty, WL)

## 
##  Moran I test under randomisation
## 
## data:  poverty  
## weights: WL    
## 
## Moran I statistic standard deviate = 6.5058, p-value = 3.864e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.442724733      -0.012345679       0.004892789

moran.plot(poverty, WL)

Local Moran

local <- localmoran(poverty,WL)

REGRESI LINEAR (OLS)

fit.ols <- lm(poverty~income+unemployment+education,data=Miss)
summary(fit.ols)

## 
## Call:
## lm(formula = poverty ~ income + unemployment + education, data = Miss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6536  -2.1611  -0.1767   2.2850   9.4530 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.409e+01  4.224e+00  10.437  < 2e-16 ***
## income       -4.425e-04  5.211e-05  -8.493 1.04e-12 ***
## unemployment  4.472e-01  1.586e-01   2.820  0.00609 ** 
## education    -2.893e-01  1.128e-01  -2.565  0.01225 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.555 on 78 degrees of freedom
## Multiple R-squared:  0.6468, Adjusted R-squared:  0.6332 
## F-statistic: 47.62 on 3 and 78 DF,  p-value: < 2.2e-16

AIC(fit.ols)

## [1] 446.6135

Pengujian Asumsi OLS

1. Normalitas Residual

shapiro.test(fit.ols$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit.ols$residuals
## W = 0.99006, p-value = 0.7868

hist(fit.ols$residuals)

2. Homoskedastisitas

bptest(fit.ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit.ols
## BP = 18.307, df = 3, p-value = 0.0003802

3. Multikolinearitas

vif(fit.ols)

##       income unemployment    education 
##     2.041353     1.308597     1.819259

4. Autokorelasi Spasial OLS

moran.test(fit.ols$residuals,listw=WL, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  fit.ols$residuals  
## weights: WL    
## 
## Moran I statistic standard deviate = 2.368, p-value = 0.01788
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.152830756      -0.012345679       0.004865573

GWR

Penentuan Bandwidth optimum

gwr.band <- gwr.sel(poverty~income+unemployment+education,coords = Coordk, Miss, adapt = TRUE)

## Adaptive q: 0.381966 CV score: 1109.623 
## Adaptive q: 0.618034 CV score: 1115.112 
## Adaptive q: 0.236068 CV score: 1116.534 
## Adaptive q: 0.4371183 CV score: 1108.453 
## Adaptive q: 0.4527045 CV score: 1109.703 
## Adaptive q: 0.4169403 CV score: 1112.953 
## Adaptive q: 0.429411 CV score: 1108.435 
## Adaptive q: 0.4309278 CV score: 1108.439 
## Adaptive q: 0.4246476 CV score: 1109.339 
## Adaptive q: 0.4275916 CV score: 1108.428 
## Adaptive q: 0.4264671 CV score: 1108.574 
## Adaptive q: 0.4282865 CV score: 1108.431 
## Adaptive q: 0.427162 CV score: 1108.426 
## Adaptive q: 0.4268966 CV score: 1108.425 
## Adaptive q: 0.4267325 CV score: 1108.465 
## Adaptive q: 0.426998 CV score: 1108.426 
## Adaptive q: 0.4268339 CV score: 1108.425 
## Adaptive q: 0.4267932 CV score: 1108.44 
## Adaptive q: 0.4268339 CV score: 1108.425

Membangun model GWR

gwr.fit <- gwr(poverty~income+unemployment+education,coords = Coordk, Miss,adapt=gwr.band,se.fit=TRUE,hatmatrix=TRUE)
gwr.fit

## Call:
## gwr(formula = poverty ~ income + unemployment + education, data = Miss, 
##     coords = Coordk, adapt = gwr.band, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.4268339 (about 35 of 82 data points)
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. 40.75790595 42.00806323 44.04947906 45.58147999 46.96809125
## income       -0.00047169 -0.00046099 -0.00043899 -0.00041259 -0.00039214
## unemployment  0.19979498  0.27891596  0.42481817  0.55455638  0.65688737
## education    -0.42805844 -0.35937280 -0.27802047 -0.18298637 -0.12847098
##               Global
## X.Intercept. 44.0917
## income       -0.0004
## unemployment  0.4472
## education    -0.2893
## Number of data points: 82 
## Effective number of parameters (residual: 2traceS - traceS'S): 8.860356 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 73.13964 
## Sigma (residual: 2traceS - traceS'S): 3.417585 
## Effective number of parameters (model: traceS): 6.908921 
## Effective degrees of freedom (model: traceS): 75.09108 
## Sigma (model: traceS): 3.372886 
## Sigma (ML): 3.227668 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 442.6204 
## AIC (GWR p. 96, eq. 4.22): 431.7835 
## Residual sum of squares: 854.263 
## Quasi-global R2: 0.693928

Uji Asumsi GWR

1. Normalitas Residual

res_gwr <- gwr.fit$SDF$gwr.e
shapiro.test(res_gwr)

## 
##  Shapiro-Wilk normality test
## 
## data:  res_gwr
## W = 0.98746, p-value = 0.6122

hist(res_gwr)

2. Homoskedastisitas

bptest(gwr.fit$lm, weights = gwr.fit$gweight)

## 
##  studentized Breusch-Pagan test
## 
## data:  gwr.fit$lm
## BP = 18.307, df = 3, p-value = 0.0003802

3. Autokorelasi Spasial GWR

gwr.morantest(gwr.fit, WL, zero.policy = TRUE)

## 
##  Leung et al. 2000 three moment approximation for Moran's I
## 
## data:  GWR residuals
## statistic = 274.63, df = 325.76, p-value = 0.01821
## sample estimates:
##          I 
## 0.08722306

Karena p-value < 0.05, maka H0 ditolak. Sehingga dapat disimpulkan bahwa terdapat autokorelasi spasial pada GWR

5. Uji Signfikansi Parameter

LMZ.F3GWR.test(gwr.fit)

## 
## Leung et al. (2000) F(3) test
## 
##              F statistic Numerator d.f. Denominator d.f.     Pr(>)    
## (Intercept)       1.1073        27.5032           75.959    0.3542    
## income            1.2157        22.0024           75.959    0.2611    
## unemployment      3.8869        31.8843           75.959 6.505e-07 ***
## education         3.8886        24.2460           75.959 3.045e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dapa dilihat, bahwa variabel income merupakan variabel satu-satunya yang tidak berpengaruh terhadap angka kemiskinan

Peta Persebaran Parameter GWR

Miss$b1 <- gwr.fit$SDF$unemployment
Miss$b2 <- gwr.fit$SDF$income 
Miss$b3 <- gwr.fit$SDF$education  

ggplot(Miss) +
  geom_sf(aes(fill = b1)) +
  scale_fill_gradient(low = "yellow", high = "red") +  
  labs(title = "Peta Sebaran Penduga Parameter b1 (unemployment") +
  theme_minimal()

ggplot(Miss) +
  geom_sf(aes(fill = b2)) +
  scale_fill_gradient(low = "yellow", high = "red") + 
  labs(title = "Peta Sebaran Penduga Parameter b2 (income)") +
  theme_minimal()

ggplot(Miss) +
  geom_sf(aes(fill = b3)) +
  scale_fill_gradient(low = "yellow", high = "red") + 
  labs(title = "Peta Sebaran Penduga Parameter b3 (education)") +
  theme_minimal()

PERBANDINGAN PENYEBARAN NILAI DUGAAN POVERTY MISSISSIPPI

gwr.pred <- gwr.fit$SDF$pred 
Miss$pred <- gwr.pred 
Miss$predols <- predict(fit.ols, data=data)

ggplot(Miss) +
  geom_sf(aes(fill = poverty)) +
  scale_fill_gradient(low = "yellow", high = "red", name = "Poverty") +  
  labs(title = "Peta Sebaran Nilai Aktual Poverty") +
  theme_minimal()

ggplot(Miss) +
  geom_sf(aes(fill = pred)) +
  scale_fill_gradient(low = "yellow", high = "red", name = "Predicted Poverty GWR") + 
  labs(title = "Peta Sebaran Nilai Dugaan Poverty GWR") +
  theme_minimal()

ggplot(Miss) +
  geom_sf(aes(fill = predols)) +
  scale_fill_gradient(low = "yellow", high = "red", name = "Predicted Poverty OLS") + 
  labs(title = "Peta Sebaran Nilai Dugaan Poverty OLS") +
  theme_minimal()

PERBANDINGAN MODEL OLS DENGAN GWR

Model dibandingkan dengan AIC. Semakin kecil AIC, maka semakin baik model

data.frame("MODEL" = c("GWR","OLS"),
           "AIC" = c(gwr.fit[["results"]][["AICh"]],AIC(fit.ols)))%>% arrange(AIC)

##   MODEL      AIC
## 1   GWR 431.7835
## 2   OLS 446.6135

Dapat dilihat bahwa model GWR lebih baik daripada OLS, karena mampu menangkap variabilitas spasial dengan baik

Aplikasi Metode Geogrephically Weighted Regression (GWR) dalam Memodelkan Tingat Kemiskinan di Negara Bagian Mississippi

Muhammad Kamil Dipinto

2024-10-16

INPUT DATA DAN PLOT PETA

PLOT PERSEBARAN TIAP VARIABEL

1. Persebaran Kemiskinan

2. Persebaran Penghasilan Rumah Tangga

3. Persebaran Pengangguran

4. Persebaran Edukasi dibawah SMA

MEMBANGUN MATRIKS PEMBOBOT SPASIAL DENGAN ROOK

Mendapatkan Koordinat Centroid

Mendapatkan Matriks Pembobot Spasial

Plot Hubungan Lokasi pada Peta

Moran’s Test

Local Moran

REGRESI LINEAR (OLS)

Pengujian Asumsi OLS

1. Normalitas Residual

2. Homoskedastisitas

3. Multikolinearitas

4. Autokorelasi Spasial OLS

GWR

Penentuan Bandwidth optimum

Membangun model GWR

Uji Asumsi GWR

1. Normalitas Residual

2. Homoskedastisitas

3. Autokorelasi Spasial GWR

5. Uji Signfikansi Parameter

Peta Persebaran Parameter GWR

PERBANDINGAN PENYEBARAN NILAI DUGAAN POVERTY MISSISSIPPI

PERBANDINGAN MODEL OLS DENGAN GWR