Penanganan Autokorelasi Spasial pada Residual OLS

Muat Package

Muat package yang dibutuhkan

library(tidyverse) 
library(sf)
library(sp)
library(MuMIn)
library(fields)
library(viridisLite)
library(gridExtra) 
library(knitr)
library(gt)
library(terra) 
library(geodata)
library(leaflet)
library(maps)
library(raster)
library(stars)
library(gtools)
library(spdep)  
library(splancs)
library(mgcv)
library(DAAG)
library(reshape2)
library(gstat)
library(oce)  
library(chirps)
library(spgwr)
library(car)
library(lmtest)
library(tmap)
library(tmaptools)
library(grid)

Muat Data

Dalam hal ini data yang akan digunakan terdiri dari Rata-rata Lama Sekolah (RLS) di Provinsi Papua pada tahun 2023 sebagai variabel dependen (Y), Pengeluaran Per Kapita (PPK) di Provinsi Papua pada tahun 2023 sebagai variabel independen pertama (X1) dan Jumlah Penduduk Miskin (JPM) di Provinsi Papua pada tahun 2023 sebagai variabel independen kedua (X2)

setwd("C:/Users/6atik/Semester 5/Spatial")
data=readxl::read_excel("datapapua.xlsx")
head(data)

## # A tibble: 6 × 4
##   `Kabupaten/kota`   rls   ppk   jpm
##   <chr>            <dbl> <dbl> <dbl>
## 1 Merauke           9    10747  24.0
## 2 Jayawijaya        5.32  8159  73.4
## 3 Jayapura          9.55 10671  14.9
## 4 Nabire            9.7   9381  34.8
## 5 Kepulauan Yapen   8.85  8091  26.0
## 6 Biak Numfor       9.8  10229  35.7

Analisis Data Eksplorasi

Analisis data eksplorasi yang dilakukan meliputi membuat ringkasan singkat dari keseluruhan data, membuat peta tematik pesebaran tiap variabel yang terkandung di dalam data berdasarkan kabupaten/kota

Melihat ringkasan data secara singkat:

summary(data)

##  Kabupaten/kota          rls              ppk             jpm       
##  Length:29          Min.   : 1.370   Min.   : 4352   Min.   : 5.45  
##  Class :character   1st Qu.: 2.710   1st Qu.: 5397   1st Qu.:19.83  
##  Mode  :character   Median : 4.900   Median : 6259   Median :26.82  
##                     Mean   : 5.769   Mean   : 7346   Mean   :31.56  
##                     3rd Qu.: 8.720   3rd Qu.: 8396   3rd Qu.:42.01  
##                     Max.   :11.670   Max.   :15272   Max.   :76.43

Membuat Peta Provinsi Papua (Sebelum Pemekaran)

Indo_Prov<-readRDS('gadm36_IDN_2_sp.rds')#memuat file 
Papua<-Indo_Prov[Indo_Prov$NAME_1 == "Papua",]#memilih Provinsi Papua saja
plot(Papua)#membuat plot Provinsi Papua

Menggabungkan data longitude dan latitude Provinsi Papua yang digunakan untuk membuat peta dengan data yang berisi variabel dependen dan independen

Papua_sf<-st_as_sf(Papua)
data=rename(data,NAME_2='Kabupaten/kota')
Papua_merged <- Papua_sf %>%
  left_join(data, by = "NAME_2") 
Papua_merged$id<-c(1:29)

Membuat Peta Sebaran RLS

#PETA KONTINU
ggplot() +
  geom_sf(data=Papua_merged, aes(fill = rls),color=NA) +
  theme_bw() +
  scale_fill_gradient(low = "#e8f6fb", high = "#0060ca") +
  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 y-axis labels   
  labs(title = "Rata-rata Lama Sekolah",
       fill = "Tahun")

#PETA DISKRIT
breaks <- c(-Inf,2.71, 4.9, 8.72,Inf) #berdasarkan kuartil

# Membuat label dari tiap interval
labels <- c("Sangat Rendah", "Rendah", "Tinggi", "Sangat Tinggi")

# Membuat kolom baru dari data yang sudah dipartisi dengan kuartil
Papua_merged$rls_Discrete <- cut(Papua_merged$rls, breaks = breaks, labels = labels, right = TRUE)

ggplot() +
  geom_sf(data=Papua_merged, aes(fill = rls_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Sangat Rendah" = "#e8f6fb", 
                               "Rendah" = "#A7C4E4",
                               "Tinggi" = "#5c9bdf",
                               "Sangat Tinggi" = "#0060ca"))+
  labs(fill = "rls_Discrete")+theme(legend.position = "right",
                             axis.text.x = element_blank(),  # Remove x-axis labels
                             axis.text.y = element_blank())+  # Remove y-axis labels   
  labs(title = "Rata rata Lama Sekolah",
       fill = "Tahun")

Membuat Peta Sebaran PPK

#PETA KONTINU
ggplot() +
  geom_sf(data=Papua_merged, aes(fill = ppk),color=NA) +
  theme_bw() +
  scale_fill_gradient(low = "#e6fff7", high = "#117555") +
  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 y-axis labels   
  labs(title = "Pengeluaran Per Kapita",
       fill = "Ribu Rupiah")

#PETA PESEBARAN DISKRIT
breaks <- c(-Inf,5397, 6259, 8396 ,Inf)#Membuat interval berdasarkan kuartil

# Menentukan label dari tiap interval
labels <- c("Sangat Rendah", "Rendah", "Tinggi", "Sangat Tinggi")

#Membuat kolom baru untuk data yang sudah di bagi
Papua_merged$ppk_Discrete <- cut(Papua_merged$ppk, breaks = breaks, labels = labels, right = TRUE)


ggplot() +
  geom_sf(data=Papua_merged, aes(fill = ppk_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Sangat Rendah" = "#e6fff7", 
                               "Rendah" = "#a5ffe3",
                               "Tinggi" = "#66ccac",
                               "Sangat Tinggi" = "#117555"))+
  labs(fill = "ppk_Discrete")+theme(legend.position = "right",
                             axis.text.x = element_blank(),  # Remove x-axis labels
                             axis.text.y = element_blank())+  # Remove y-axis labels   
  labs(title = "Pengeluaran Per Kapita",
       fill = "Ribu Rupiah")

Membuat Peta Sebaran JPM

#Jumlah Penduduk Miskin
#Maps KONTINU
ggplot() +
  geom_sf(data=Papua_merged, aes(fill = jpm),color=NA) +
  theme_bw() +
  scale_fill_gradient(low = "#f5f7bb", high = "#72770a") +
  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 y-axis labels   
  labs(title = "Jumlah Penduduk Miskin",
       fill = "Ribu Jiwa")

#PETA DISKRIT
breaks <- c(-Inf,19.83, 26.82, 42.01 ,Inf) #membagi data ke dalam interval berdasarkan kuartil

# Menentukan label dari tiap interval
labels <- c("Sangat Rendah", "Rendah", "Tinggi", "Sangat Tinggi")

# Membuat kolom baru untuk data yang sudah diintervalkan
Papua_merged$jpm_Discrete <- cut(Papua_merged$jpm, breaks = breaks, labels = labels, right = TRUE)


ggplot() +
  geom_sf(data = Papua_merged, aes(fill = jpm_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Sangat Rendah" = "#f5f7bb", 
                               "Rendah" = "#d6dc5d",
                               "Tinggi" = "#adb411",
                               "Sangat Tinggi" = "#72770a"))+
  labs(fill = "ppk_Discrete")+theme(legend.position = "right",
                             axis.text.x = element_blank(),  # Remove x-axis labels
                             axis.text.y = element_blank())+  # Remove y-axis labels   
  labs(title = "Jumlah Penduduk Miskin",
       fill = "Ribu Jiwa")

Membuat Model Ordinary Least Square (OLS)

Model OLS nantinya akan dibandingkan dengan model GWR agar dapat menentukan model yang lebih tepat

model_ols <- lm(rls~ppk+jpm, data = Papua_merged)
summary(model_ols)

## 
## Call:
## lm(formula = rls ~ ppk + jpm, data = Papua_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9848 -0.6989 -0.3179  1.1165  2.5734 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.1177543  1.0584450  -0.111  0.91227    
## ppk          0.0010005  0.0001136   8.809 2.78e-09 ***
## jpm         -0.0463462  0.0147302  -3.146  0.00411 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 26 degrees of freedom
## Multiple R-squared:  0.7925, Adjusted R-squared:  0.7766 
## F-statistic: 49.66 on 2 and 26 DF,  p-value: 1.32e-09

Mengecek Asumsi Model OLS

Model OLS yang digunakan harus memenuhi beberapa asumsi, yaitu:

Homoskedastisitas

bptest(model_ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_ols
## BP = 1.4842, df = 2, p-value = 0.4761

Model memnuhi asumsi homoskedastisitas

Multikolinearitas

Gunakan VIF

vif(model_ols)

##      ppk      jpm 
## 1.027896 1.027896

Model memenuhi asumsi non multikolinearitas

Normalitas Residual

shapiro.test(model_ols$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_ols$residuals
## W = 0.97199, p-value = 0.6146

Model Memenuhi asumsi normalitas residual

Autokorelasi Temporal

dwtest(model_ols)

## 
##  Durbin-Watson test
## 
## data:  model_ols
## DW = 1.978, p-value = 0.483
## alternative hypothesis: true autocorrelation is greater than 0

Model memenuhi asumsi non autokorelasi temporal

Autokorelasi Spasial

CoordK=coordinates(Papua)
# Membuat matriks bobot spasial berdasarkan tetangga terdekat
nb = knn2nb(knearneigh(CoordK, k=3))
listw = nb2listw(nb, style="W")

# Uji Moran's I untuk residual spasial
moran.test(model_ols$residuals, listw)

## 
##  Moran I test under randomisation
## 
## data:  model_ols$residuals  
## weights: listw    
## 
## Moran I statistic standard deviate = 2.5939, p-value = 0.004745
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.29715473       -0.03571429        0.01646839

Uji autokorelasi spasial mengindikasikan terdapat nya autokorelasi spasial pada residual model OLS, hal ini menyebabkan residual OLS tidak lagi iid sehingga salah satu asumsi dalam OLS terlanggar. Untuk menangani pelanggaran asumsi OLS dibuatlah model GWR yang dapat menangani pola spasial pada data

Membuat Model Geographically Weighted Regression (GWR)

Menentukan Kernel dan Bandwith yang tepat

Jenis Kernel di tentukan menggunakan nilai CV, semakin kecil nilai CV semakin bagus kernel

Kernel Bisquare Fixed memiliki nilai CV terkecil sebesar 60.79 dengan bandwith 3.28

bandwith_bifix = gwr.sel(rls~ppk+jpm, 
                   data=data,
                   coords = coordinates(Papua),
                   gweight=gwr.bisquare,
                   adapt=F)

## Bandwidth: 3.361579 CV score: 60.92438 
## Bandwidth: 5.433718 CV score: 69.13346 
## Bandwidth: 2.080926 CV score: NA

## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value

## Bandwidth: 4.153066 CV score: 67.91251 
## Bandwidth: 2.872413 CV score: NA

## Warning in optimize(gwr.cv.f, lower = beta1, upper = beta2, maximum = FALSE, :
## NA/Inf replaced by maximum positive value

## Bandwidth: 3.6639 CV score: 66.71158 
## Bandwidth: 3.174734 CV score: 61.1446 
## Bandwidth: 3.282342 CV score: 60.79572 
## Bandwidth: 3.290786 CV score: 60.79681 
## Bandwidth: 3.283506 CV score: 60.79567 
## Bandwidth: 3.283925 CV score: 60.79566 
## Bandwidth: 3.283884 CV score: 60.79566 
## Bandwidth: 3.283965 CV score: 60.79566 
## Bandwidth: 3.283925 CV score: 60.79566

Kernel Bisquare Adaptive

Kernel Bisquare Adaptive memiliki nilai CV terkecil sebesar 65.28 dengan bandwith 0.276

bandwith_biadapt = gwr.sel(rls~ppk+jpm, 
                   data=data,
                   coords = coordinates(Papua),
                   gweight=gwr.bisquare,
                   adapt=T)

## Adaptive q: 0.381966 CV score: 90.04317 
## Adaptive q: 0.618034 CV score: 68.25351 
## Adaptive q: 0.763932 CV score: 65.8354 
## Adaptive q: 0.7327817 CV score: 65.3597 
## Adaptive q: 0.7208436 CV score: 65.41452 
## Adaptive q: 0.7317933 CV score: 65.35017 
## Adaptive q: 0.7285795 CV score: 65.32053 
## Adaptive q: 0.7256247 CV score: 65.29432 
## Adaptive q: 0.7237985 CV score: 65.29456 
## Adaptive q: 0.7247469 CV score: 65.28661 
## Adaptive q: 0.7247062 CV score: 65.28626 
## Adaptive q: 0.7243595 CV score: 65.28322 
## Adaptive q: 0.7241452 CV score: 65.28134 
## Adaptive q: 0.7240128 CV score: 65.28616 
## Adaptive q: 0.7242187 CV score: 65.28198 
## Adaptive q: 0.7240946 CV score: 65.28297 
## Adaptive q: 0.7241452 CV score: 65.28134

Kernel Gaussian Fixed

bandwith_gausadapt = gwr.sel(rls~ppk+jpm, 
                   data=data,
                   coords = coordinates(Papua),
                   gweight=gwr.Gauss,
                   adapt=F)

## Bandwidth: 3.361579 CV score: 74.40142 
## Bandwidth: 5.433718 CV score: 76.10512 
## Bandwidth: 2.080926 CV score: 69.50965 
## Bandwidth: 1.289439 CV score: 64.09978 
## Bandwidth: 0.8002734 CV score: 68.13802 
## Bandwidth: 1.395153 CV score: 64.40676 
## Bandwidth: 1.264895 CV score: 64.07827 
## Bandwidth: 1.249018 CV score: 64.07554 
## Bandwidth: 1.252021 CV score: 64.07537 
## Bandwidth: 1.252137 CV score: 64.07537 
## Bandwidth: 1.252096 CV score: 64.07537 
## Bandwidth: 1.252177 CV score: 64.07537 
## Bandwidth: 1.252137 CV score: 64.07537

Kernel Gaussian Adaptive

Kernel Gaussian Adaptive memiliki nilai CV terkecil sebesar 69.61 dengan bandwith 0.222

bandwith_gausadapt = gwr.sel(rls~ppk+jpm, 
                   data=data,
                   coords = coordinates(Papua),
                   gweight=gwr.Gauss,
                   adapt=T)

## Adaptive q: 0.381966 CV score: 70.42437 
## Adaptive q: 0.618034 CV score: 73.38411 
## Adaptive q: 0.236068 CV score: 69.71288 
## Adaptive q: 0.1874461 CV score: 70.56686 
## Adaptive q: 0.2878809 CV score: 69.76953 
## Adaptive q: 0.2483791 CV score: 69.72412 
## Adaptive q: 0.2174961 CV score: 69.62904 
## Adaptive q: 0.206018 CV score: 69.92125 
## Adaptive q: 0.2245899 CV score: 69.61034 
## Adaptive q: 0.2231584 CV score: 69.60746 
## Adaptive q: 0.2226492 CV score: 69.60715 
## Adaptive q: 0.2224754 CV score: 69.60713 
## Adaptive q: 0.2225161 CV score: 69.60713 
## Adaptive q: 0.2225568 CV score: 69.60713 
## Adaptive q: 0.2225161 CV score: 69.60713

Kernel dengan merupakan nilai CV terkecil merupakan kernel Bisquare Fixed dengan bandwtih 3.28

Sehingga model GWR yang akan digunakan adalah berikut:

rtg.model <- gwr(rls~ppk+jpm,
                 data=Papua_merged, 
                 coords = coordinates(Papua),
                 gweight=gwr.bisquare,
                 bandwidth = bandwith_bifix,
                 hatmatrix=TRUE)

rtg.model

## Call:
## gwr(formula = rls ~ ppk + jpm, data = Papua_merged, coords = coordinates(Papua), 
##     bandwidth = bandwith_bifix, gweight = gwr.bisquare, hatmatrix = TRUE)
## Kernel function: gwr.bisquare 
## Fixed bandwidth: 3.283925 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -2.96812856 -1.93927223 -1.08682162 -0.60104409  4.90850627
## ppk           0.00065678  0.00096750  0.00109382  0.00124482  0.00155047
## jpm          -0.13044830 -0.04527098 -0.03627370 -0.03008081 -0.02238339
##               Global
## X.Intercept. -0.1178
## ppk           0.0010
## jpm          -0.0463
## Number of data points: 29 
## Effective number of parameters (residual: 2traceS - traceS'S): 11.30889 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 17.69111 
## Sigma (residual: 2traceS - traceS'S): 1.298578 
## Effective number of parameters (model: traceS): 9.449376 
## Effective degrees of freedom (model: traceS): 19.55062 
## Sigma (model: traceS): 1.23528 
## Sigma (ML): 1.014254 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 117.6516 
## AIC (GWR p. 96, eq. 4.22): 92.56869 
## Residual sum of squares: 29.83261 
## Quasi-global R2: 0.8977016

Uji Asumsi GWR

Model GWR yang digunakan harus memenuhi beberapa asumsi, yaitu

Non Autokorelasi Spasial

Uji Autokorelasi Spatial gilakuakn menggunakan uji Moran’s, namun sebelumnya perlu ditentukan matriks bobot yang sesuai

Penentuan Matriks Bobot

Matriks Bobot Rook

CoordK=coordinates(Papua)
WR <- poly2nb(Papua_merged, queen=FALSE) ; WR

## Neighbour list object:
## Number of regions: 29 
## Number of nonzero links: 120 
## Percentage nonzero weights: 14.26873 
## Average number of links: 4.137931 
## 1 region with no links:
## 10
## 3 disjoint connected subgraphs

WBR <- nb2mat(WR, style='B', zero.policy= TRUE) #Dalam bentuk matrix biner "B"
WBR[c(1:5),c(1:5)]

##   [,1] [,2] [,3] [,4] [,5]
## 1    0    0    0    0    0
## 2    0    0    0    0    0
## 3    0    0    0    0    0
## 4    0    0    0    0    1
## 5    0    0    0    1    0

WLR<-nb2listw(WR, zero.policy = T) ; # Moran's I (List Neighbours)
plot(Papua, axes=T, col="gray90")
plot(WR,coordinates(Papua),col="red",add=TRUE)
text(CoordK[,1], CoordK[,2], row.names(Papua), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue") #hubungan antar lokasi

ujimoranr=moran.test(Papua_merged$rls,WLR);print(ujimoranr)

## 
##  Moran I test under randomisation
## 
## data:  Papua_merged$rls  
## weights: WLR  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.7561, p-value = 0.03953
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.20471122       -0.03703704        0.01895064

plotmoranr=moran.plot(Papua_merged$rls,WLR); plotmoranr

##        x       wx is_inf labels       dfb.1_        dfb.x        dffit
## 1   4.83 5.305000  FALSE      1  0.004122162 -0.001647208  0.005799811
## 2   9.80 8.350000  FALSE      2 -0.148645254  0.298855188  0.380258947
## 3   8.68 5.265000  FALSE      3  0.015972802 -0.045433827 -0.067189093
## 4   1.83 5.920000  FALSE      4  0.166886539 -0.132020320  0.169481050
## 5   4.13 6.290000  FALSE      5  0.098144935 -0.054279921  0.118200489
## 6   2.06 5.812500  FALSE      6  0.138701177 -0.107675509  0.141660481
## 7   9.55 6.260000  FALSE      7 -0.018078293  0.038275263  0.049956616
## 8   5.32 2.624000  FALSE      8 -0.151880372  0.035519657 -0.253192423
## 9   7.93 7.976667  FALSE      9 -0.024782611  0.146293043  0.259807942
## 10  8.85 0.000000   TRUE     10  0.241090238 -0.633579858 -0.909262607
## 11 11.67 8.740000  FALSE     11 -0.306811223  0.496876429  0.564083739
## 12  2.21 2.893333  FALSE     12 -0.277918880  0.212883952 -0.285118796
## 13  4.90 4.628333  FALSE     13 -0.040860996  0.015516085 -0.058688617
## 14  2.67 4.470000  FALSE     14 -0.057866822  0.042263566 -0.060466056
## 15  6.15 7.503333  FALSE     15  0.071634533  0.022842241  0.191640346
## 16  9.00 7.415000  FALSE     16 -0.061183078  0.151615280  0.212448684
## 17 10.17 2.886667  FALSE     17  0.233053868 -0.440351992 -0.542769800
## 18  9.70 4.592500  FALSE     18  0.078563878 -0.161077009 -0.206961009
## 19  1.70 5.146000  FALSE     19  0.054250194 -0.043340863  0.054946767
## 20  3.46 4.876667  FALSE     20 -0.017478959  0.011366297 -0.019308138
## 21  2.71 7.340000  FALSE     21  0.305173048 -0.221841883  0.319516077
## 22  1.37 5.068571  FALSE     22  0.053007522 -0.043322613  0.053404382
## 23  3.96 2.705000  FALSE     23 -0.221845528  0.128936687 -0.260178687
## 24  8.71 5.115000  FALSE     24  0.021167036 -0.059305444 -0.087218451
## 25  8.35 9.800000  FALSE     25 -0.093720678  0.329472814  0.521974371
## 26  2.34 3.435000  FALSE     26 -0.196916834  0.148980419 -0.202913733
## 27  8.72 4.507500  FALSE     27  0.040512460 -0.112950040 -0.165809254
## 28  3.20 5.161429  FALSE     28  0.016749687 -0.011387730  0.018089367
## 29  3.34 5.890000  FALSE     29  0.091666633 -0.060903403  0.100147788
##        cov.r       cook.d        hat
## 1  1.1203533 0.0000174652 0.03750825
## 2  1.0627271 0.0710909560 0.09019317
## 3  1.1456947 0.0023380185 0.06353432
## 4  1.1553551 0.0147448334 0.08769580
## 5  1.1016176 0.0071700395 0.04369786
## 6  1.1541772 0.0103301022 0.08166342
## 7  1.1741731 0.0012944614 0.08349671
## 8  0.9805966 0.0311774708 0.03517502
## 9  1.0324736 0.0334167124 0.05049169
## 10 0.5554073 0.2975699073 0.06702696
## 11 1.1188697 0.1547974049 0.15387695
## 12 1.0875751 0.0407037584 0.07792471
## 13 1.1122553 0.0017822820 0.03707412
## 14 1.1518810 0.0018946971 0.06742156
## 15 1.0352354 0.0183540397 0.03497972
## 16 1.1084152 0.0229090451 0.07027329
## 17 0.9894919 0.1389356670 0.10089057
## 18 1.1421624 0.0218643312 0.08746315
## 19 1.1839689 0.0015658237 0.09126588
## 20 1.1378938 0.0001935216 0.05276975
## 21 1.0378990 0.0502426194 0.06657682
## 22 1.1970136 0.0014794138 0.10084894
## 23 1.0165392 0.0333361996 0.04570821
## 24 1.1425307 0.0039330276 0.06413626
## 25 0.8323598 0.1206718522 0.05732029
## 26 1.1212505 0.0209681365 0.07480932
## 27 1.1179095 0.0140588662 0.06433828
## 28 1.1432565 0.0001698701 0.05711934
## 29 1.1257740 0.0051731928 0.05471964

locmR<-localmoran(Papua_merged$rls,WLR) ; locmR

##             Ii          E.Ii      Var.Ii         Z.Ii Pr(z != E(Ii))
## 1   0.04337037 -0.0031335442 0.020130673  0.327763320    0.743090623
## 2   1.03440543 -0.0577000699 1.576752384  0.869726786    0.384449731
## 3  -0.14597196 -0.0300891156 0.188073125 -0.267211880    0.789306029
## 4  -0.05903088 -0.0551135024 0.466112384 -0.005737868    0.995421869
## 5  -0.08488203 -0.0095442098 0.060920093 -0.305233921    0.760187999
## 6  -0.01593120 -0.0488656814 0.299523772  0.060177640    0.952014154
## 7   0.18448182 -0.0507644489 0.189775147  0.540011590    0.589189039
## 8   0.14053535 -0.0007169856 0.003539900  2.374106059    0.017591497
## 9   0.47428620 -0.0165806825 0.145946648  1.284892126    0.198829983
## 10  0.00000000  0.0000000000 0.000000000          NaN            NaN
## 11  1.74315380 -0.1236582697 1.513123040  1.517621526    0.129109837
## 12  1.01795110 -0.0449934538 0.169223762  2.583925211    0.009768297
## 13  0.09863417 -0.0026839061 0.010541582  0.986811107    0.323735249
## 14  0.40045490 -0.0341151865 0.129771329  1.206342010    0.227685646
## 15  0.06564499 -0.0005147066 0.004604571  0.974987195    0.329566612
## 16  0.52871159 -0.0370687636 0.498403360  0.801414842    0.422891526
## 17 -1.26150029 -0.0687795140 0.252241935 -2.374816950    0.017557655
## 18 -0.45999264 -0.0548725483 0.334218889 -0.700758358    0.483453830
## 19  0.25223236 -0.0588110919 0.273481598  0.594780886    0.551989946
## 20  0.20499161 -0.0189400979 0.073178484  0.827797071    0.407785434
## 21 -0.47784777 -0.0332402812 0.207094574 -0.976995174    0.328571541
## 22  0.30656071 -0.0687364056 0.206259961  0.826357085    0.408601560
## 23  0.55134249 -0.0116263596 0.074054318  2.068755521    0.038569034
## 24 -0.19134126 -0.0307125591 0.191846586 -0.366729722    0.713820634
## 25  1.03440543 -0.0236531606 0.669716968  1.292897411    0.196046543
## 26  0.79605126 -0.0417667995 0.257921708  1.649703020    0.099003692
## 27 -0.37024847 -0.0309217943 0.193111882 -0.772171020    0.440013143
## 28  0.15531426 -0.0234450276 0.073773932  0.658138282    0.510449279
## 29 -0.02915605 -0.0209596303 0.101385602 -0.025741632    0.979463417
## attr(,"call")
## localmoran(x = Papua_merged$rls, listw = WLR)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1    Low-Low  Low-High   Low-Low
## 2  High-High High-High High-High
## 3   High-Low High-High  High-Low
## 4   Low-High  Low-High  Low-High
## 5   Low-High  Low-High  Low-High
## 6   Low-High  Low-High  Low-High
## 7  High-High High-High High-High
## 8    Low-Low  High-Low   Low-Low
## 9  High-High High-High High-High
## 10  High-Low  High-Low  High-Low
## 11 High-High High-High High-High
## 12   Low-Low   Low-Low   Low-Low
## 13   Low-Low   Low-Low   Low-Low
## 14   Low-Low   Low-Low   Low-Low
## 15 High-High High-High High-High
## 16 High-High High-High High-High
## 17  High-Low  High-Low  High-Low
## 18  High-Low  High-Low  High-Low
## 19   Low-Low   Low-Low   Low-Low
## 20   Low-Low   Low-Low   Low-Low
## 21  Low-High  Low-High  Low-High
## 22   Low-Low   Low-Low   Low-Low
## 23   Low-Low   Low-Low   Low-Low
## 24  High-Low  High-Low  High-Low
## 25 High-High High-High High-High
## 26   Low-Low   Low-Low   Low-Low
## 27  High-Low  High-Low  High-Low
## 28   Low-Low   Low-Low   Low-Low
## 29  Low-High  Low-High  Low-High

Matriks bobot Rook menghasilkan p-value sebesar 0.03953

Matriks bobot Queen

WQ <- poly2nb(Papua_merged, queen=TRUE) #mendapatkan W
WBQ <- nb2mat(WQ, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"  
WBQ[1:5,1:5]

##   [,1] [,2] [,3] [,4] [,5]
## 1    0    0    1    0    0
## 2    0    0    0    0    0
## 3    1    0    0    0    0
## 4    0    0    0    0    1
## 5    0    0    0    1    0

WLQ<-nb2listw(WQ, zero.policy = T) ; WLQ # Moran's I (List Neighbours)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 29 
## Number of nonzero links: 128 
## Percentage nonzero weights: 15.21998 
## Average number of links: 4.413793 
## 1 region with no links:
## 10
## 3 disjoint connected subgraphs
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 28 784 28 15.52919 114.9378

plot(Papua, axes=T, col="gray90")
plot(WQ,coordinates(Papua),col="red",add=TRUE)
text(CoordK[,1], CoordK[,2], row.names(Papua), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue") #hubungan antar lokasi

ujimoranq=moran.test(Papua_merged$rls,WLQ);print(ujimoranq)

## 
##  Moran I test under randomisation
## 
## data:  Papua_merged$rls  
## weights: WLQ  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.7693, p-value = 0.03842
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.20016313       -0.03703704        0.01797401

plotq=moran.plot(Papua_merged$rls,WLQ);plotq

##        x       wx is_inf labels      dfb.1_       dfb.x       dffit     cov.r
## 1   4.83 5.980000  FALSE      1  0.05142973 -0.02055122  0.07236074 1.1089366
## 2   9.80 8.350000  FALSE      2 -0.15452999  0.31068660  0.39531306 1.0536689
## 3   8.68 5.178000  FALSE      3  0.01810513 -0.05149913 -0.07615867 1.1440305
## 4   1.83 5.920000  FALSE      4  0.17044193 -0.13483291  0.17309172 1.1542253
## 5   4.13 6.290000  FALSE      5  0.10116838 -0.05595206  0.12184176 1.1000172
## 6   2.06 5.812500  FALSE      6  0.14177150 -0.11005904  0.14479631 1.1532886
## 7   9.55 6.260000  FALSE      7 -0.02031984  0.04302106  0.05615080 1.1735225
## 8   5.32 2.624000  FALSE      8 -0.15466891  0.03617180 -0.25784105 0.9760090
## 9   7.93 7.976667  FALSE      9 -0.02567487  0.15156011  0.26916195 1.0254636
## 10  8.85 0.000000   TRUE     10  0.24716186 -0.64953595 -0.93216150 0.5384247
## 11 11.67 8.740000  FALSE     11 -0.32000107  0.51823720  0.58833376 1.1065800
## 12  2.21 3.932857  FALSE     12 -0.13261923  0.10158542 -0.13605494 1.1500789
## 13  4.90 4.628333  FALSE     13 -0.04096072  0.01555395 -0.05883186 1.1122180
## 14  2.67 4.470000  FALSE     14 -0.05889742  0.04301627 -0.06154295 1.1517203
## 15  6.15 6.427500  FALSE     15  0.03530688  0.01125837  0.09445476 1.0966309
## 16  9.00 7.415000  FALSE     16 -0.06386058  0.15825029  0.22174589 1.1039772
## 17 10.17 2.790000  FALSE     17  0.24398799 -0.46101186 -0.56823478 0.9722889
## 18  9.70 4.592500  FALSE     18  0.07830272 -0.16054157 -0.20627305 1.1424186
## 19  1.70 4.516667  FALSE     19 -0.04391359  0.03508288 -0.04447744 1.1849123
## 20  3.46 4.876667  FALSE     20 -0.01726676  0.01122831 -0.01907373 1.1379080
## 21  2.71 7.340000  FALSE     21  0.31277939 -0.22737122  0.32747992 1.0324539
## 22  1.37 4.647500  FALSE     22 -0.01718786  0.01404749 -0.01731654 1.1991096
## 23  3.96 2.705000  FALSE     23 -0.22633608  0.13154660 -0.26544518 1.0122664
## 24  8.71 5.115000  FALSE     24  0.02070807 -0.05801952 -0.08532729 1.1429475
## 25  8.35 9.800000  FALSE     25 -0.09700330  0.34101280  0.54025684 0.8151870
## 26  2.34 3.435000  FALSE     26 -0.20135697  0.15233967 -0.20748909 1.1192886
## 27  8.72 4.507500  FALSE     27  0.04047690 -0.11285089 -0.16566371 1.1179690
## 28  3.20 5.285000  FALSE     28  0.03155146 -0.02145112  0.03407502 1.1420471
## 29  3.34 5.890000  FALSE     29  0.09432289 -0.06266822  0.10304980 1.1248982
##          cook.d        hat
## 1  0.0027047542 0.03750825
## 2  0.0765030973 0.09019317
## 3  0.0030017425 0.06353432
## 4  0.0153722593 0.08769580
## 5  0.0076130669 0.04369786
## 6  0.0107883479 0.08166342
## 7  0.0016349127 0.08349671
## 8  0.0322570994 0.03517502
## 9  0.0357443033 0.05049169
## 10 0.3079281311 0.06702696
## 11 0.1674656447 0.15387695
## 12 0.0095311569 0.07792471
## 13 0.0017909625 0.03707412
## 14 0.0019626497 0.06742156
## 15 0.0045889797 0.03497972
## 16 0.0249080001 0.07027329
## 17 0.1509487670 0.10089057
## 18 0.0217216494 0.08746315
## 19 0.0010263869 0.09126588
## 20 0.0001888525 0.05276975
## 21 0.0526397691 0.06657682
## 22 0.0001556818 0.10084894
## 23 0.0346264276 0.04570821
## 24 0.0037650032 0.06413626
## 25 0.1279326092 0.05732029
## 26 0.0219051993 0.07480932
## 27 0.0140345686 0.06433828
## 28 0.0006024385 0.05711934
## 29 0.0054752161 0.05471964

locmQ<-localmoran(Papua_merged$rls,WLQ) ; locmQ

##             Ii          E.Ii      Var.Ii         Z.Ii Pr(z != E(Ii))
## 1  -0.01968013 -0.0031335442 0.015433516 -0.133191332    0.894042066
## 2   1.03440543 -0.0577000699 1.576752384  0.869726786    0.384449731
## 3  -0.17115399 -0.0300891156 0.144189396 -0.371494362    0.710269355
## 4  -0.05903088 -0.0551135024 0.466112384 -0.005737868    0.995421869
## 5  -0.08488203 -0.0095442098 0.060920093 -0.305233921    0.760187999
## 6  -0.01593120 -0.0488656814 0.299523772  0.060177640    0.952014154
## 7   0.18448182 -0.0507644489 0.189775147  0.540011590    0.589189039
## 8   0.14053535 -0.0007169856 0.003539900  2.374106059    0.017591497
## 9   0.47428620 -0.0165806825 0.145946648  1.284892126    0.198829983
## 10  0.00000000  0.0000000000 0.000000000          NaN            NaN
## 11  1.74315380 -0.1236582697 1.513123040  1.517621526    0.129109837
## 12  0.65001200 -0.0449934538 0.138455805  1.867809698    0.061788601
## 13  0.09863417 -0.0026839061 0.010541582  0.986811107    0.323735249
## 14  0.40045490 -0.0341151865 0.129771329  1.206342010    0.227685646
## 15  0.02491712 -0.0005147066 0.003315291  0.441689076    0.658714216
## 16  0.52871159 -0.0370687636 0.498403360  0.801414842    0.422891526
## 17 -1.30380348 -0.0687795140 0.206379765 -2.718578196    0.006556316
## 18 -0.45999264 -0.0548725483 0.334218889 -0.700758358    0.483453830
## 19  0.50690201 -0.0588110919 0.217992578  1.211644968    0.225648313
## 20  0.20499161 -0.0189400979 0.073178484  0.827797071    0.407785434
## 21 -0.47784777 -0.0332402812 0.207094574 -0.976995174    0.328571541
## 22  0.49077190 -0.0687364056 0.171883301  1.349552227    0.177159657
## 23  0.55134249 -0.0116263596 0.074054318  2.068755521    0.038569034
## 24 -0.19134126 -0.0307125591 0.191846586 -0.366729722    0.713820634
## 25  1.03440543 -0.0236531606 0.669716968  1.292897411    0.196046543
## 26  0.79605126 -0.0417667995 0.257921708  1.649703020    0.099003692
## 27 -0.37024847 -0.0309217943 0.193111882 -0.772171020    0.440013143
## 28  0.12374167 -0.0234450276 0.061478277  0.593618883    0.552767048
## 29 -0.02915605 -0.0209596303 0.101385602 -0.025741632    0.979463417
## attr(,"call")
## localmoran(x = Papua_merged$rls, listw = WLQ)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1   Low-High  Low-High  Low-High
## 2  High-High High-High High-High
## 3   High-Low  High-Low  High-Low
## 4   Low-High  Low-High  Low-High
## 5   Low-High  Low-High  Low-High
## 6   Low-High  Low-High  Low-High
## 7  High-High High-High High-High
## 8    Low-Low  High-Low   Low-Low
## 9  High-High High-High High-High
## 10  High-Low  High-Low  High-Low
## 11 High-High High-High High-High
## 12   Low-Low   Low-Low   Low-Low
## 13   Low-Low   Low-Low   Low-Low
## 14   Low-Low   Low-Low   Low-Low
## 15 High-High High-High High-High
## 16 High-High High-High High-High
## 17  High-Low  High-Low  High-Low
## 18  High-Low  High-Low  High-Low
## 19   Low-Low   Low-Low   Low-Low
## 20   Low-Low   Low-Low   Low-Low
## 21  Low-High  Low-High  Low-High
## 22   Low-Low   Low-Low   Low-Low
## 23   Low-Low   Low-Low   Low-Low
## 24  High-Low  High-Low  High-Low
## 25 High-High High-High High-High
## 26   Low-Low   Low-Low   Low-Low
## 27  High-Low  High-Low  High-Low
## 28   Low-Low  Low-High   Low-Low
## 29  Low-High  Low-High  Low-High

Matriks bobt Queen menghasilkan p-value sebesar 0.03842

Matriks Bobot K-Nearest Neighbors

WK <- knn2nb(knearneigh(CoordK, k = 3)) ; WK

## Neighbour list object:
## Number of regions: 29 
## Number of nonzero links: 87 
## Percentage nonzero weights: 10.34483 
## Average number of links: 3 
## Non-symmetric neighbours list

WBK <- nb2mat(WK, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WBK[1:5,1:5]

##   [,1] [,2] [,3] [,4] [,5]
## 1    0    0    0    0    0
## 2    0    0    0    0    0
## 3    0    0    0    0    0
## 4    0    0    0    0    0
## 5    0    0    0    1    0

WLK<-nb2listw(WK,zero.policy = T) ; WLK # Moran's I (List Neighbours)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 29 
## Number of nonzero links: 87 
## Percentage nonzero weights: 10.34483 
## Average number of links: 3 
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 29 841 29 15.88889 124.6667

plot(Papua, axes=T, col="gray90")
plot(WK,coordinates(Papua),col="red",add=TRUE)
text(CoordK[,1], CoordK[,2], row.names(Papua), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue") #hubungan antar lokasi

ujimorank=moran.test(Papua_merged$rls,WLK);ujimorank

## 
##  Moran I test under randomisation
## 
## data:  Papua_merged$rls  
## weights: WLK    
## 
## Moran I statistic standard deviate = 2.9593, p-value = 0.001542
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.34924544       -0.03571429        0.01692163

plotmorank=moran.plot(Papua_merged$rls,WLK);plotmorank

##        x       wx is_inf labels       dfb.1_        dfb.x        dffit
## 1   4.83 3.406667  FALSE      1 -0.079498665  0.031767504 -0.111853238
## 2   9.80 8.640000  FALSE      2 -0.164113947  0.329955400  0.419830400
## 3   8.68 4.020000  FALSE      3  0.060273776 -0.171445703 -0.253539754
## 4   1.83 5.230000  FALSE      4  0.296449308 -0.234514616  0.301058073
## 5   4.13 4.996667  FALSE      5  0.069661858 -0.038527104  0.083897000
## 6   2.06 2.220000  FALSE      6 -0.204279993  0.158585189 -0.208638475
## 7   9.55 7.646667  FALSE      7 -0.088802380  0.188011918  0.245391890
## 8   5.32 2.740000  FALSE      8 -0.116471865  0.027238811 -0.194164612
## 9   7.93 7.976667  FALSE      9 -0.028540057  0.168473442  0.299199041
## 10  8.85 8.956667  FALSE     10 -0.121100892  0.318250488  0.456727380
## 11 11.67 6.940000   TRUE     11 -0.007803502  0.012637662  0.014347027
## 12  2.21 3.120000  FALSE     12 -0.068588135  0.052538040 -0.070365016
## 13  4.90 5.003333  FALSE     13  0.034706513 -0.013179052  0.049848939
## 14  2.67 3.666667  FALSE     14 -0.009429628  0.006887016 -0.009853184
## 15  6.15 7.503333  FALSE     15  0.097589183  0.031118449  0.261075544
## 16  9.00 6.553333  FALSE     16 -0.025438879  0.063039045  0.088332536
## 17 10.17 2.450000  FALSE     17  0.331388437 -0.626153771 -0.771785669
## 18  9.70 3.140000  FALSE     18  0.201932642 -0.414016047 -0.531951638
## 19  1.70 3.290000  FALSE     19 -0.017957578  0.014346436 -0.018188154
## 20  3.46 4.686667  FALSE     20  0.075687811 -0.049218616  0.083608568
## 21  2.71 4.823333  FALSE     21  0.145350332 -0.105660678  0.152181747
## 22  1.37 3.160000  FALSE     22 -0.021815702  0.017829795 -0.021979033
## 23  3.96 1.973333  FALSE     23 -0.221279036  0.128607442 -0.259514310
## 24  8.71 5.706667  FALSE     24  0.004649326 -0.013026403 -0.019157477
## 25  8.35 9.123333  FALSE     25 -0.084476184  0.296974016  0.470487454
## 26  2.34 2.946667  FALSE     26 -0.098838802  0.074777995 -0.101848835
## 27  8.72 2.296667  FALSE     27  0.128951319 -0.359520418 -0.527771504
## 28  3.20 3.790000  FALSE     28 -0.015697019  0.010672046 -0.016952504
## 29  3.34 3.730000  FALSE     29 -0.027729502  0.018423509 -0.030295084
##        cov.r       cook.d        hat
## 1  1.0932619 6.416937e-03 0.03750825
## 2  1.0384306 8.566060e-02 0.09019317
## 3  1.0720103 3.220383e-02 0.06353432
## 4  1.1007773 4.541404e-02 0.08769580
## 5  1.1144354 3.633188e-03 0.04369786
## 6  1.1312999 2.218445e-02 0.08166342
## 7  1.1190293 3.049146e-02 0.08349671
## 8  1.0338458 1.882619e-02 0.03517502
## 9  1.0018171 4.365500e-02 0.05049169
## 10 0.9353091 9.743079e-02 0.06702696
## 11 1.2744110 1.068723e-04 0.15387695
## 12 1.1642853 2.565054e-03 0.07792471
## 13 1.1143831 1.287050e-03 0.03707412
## 14 1.1562469 5.040704e-05 0.06742156
## 15 0.9718355 3.300403e-02 0.03497972
## 16 1.1507572 4.035347e-03 0.07027329
## 17 0.8271724 2.568430e-01 0.10089057
## 18 0.9530356 1.319453e-01 0.08746315
## 19 1.1864077 1.717444e-04 0.09126588
## 20 1.1275698 3.612194e-03 0.05276975
## 21 1.1269954 1.187669e-02 0.06657682
## 22 1.1989589 2.507874e-04 0.10084894
## 23 1.0170740 3.317489e-02 0.04570821
## 24 1.1518325 1.905231e-04 0.06413626
## 25 0.8802268 1.008197e-01 0.05732029
## 26 1.1541822 5.359632e-03 0.07480932
## 27 0.8627728 1.251323e-01 0.06433828
## 28 1.1433143 1.491932e-04 0.05711934
## 29 1.1394380 4.762555e-04 0.05471964

locmK<-localmoran(Papua_merged$rls,WLK) ; locmK

##             Ii          E.Ii      Var.Ii        Z.Ii Pr(z != E(Ii))
## 1   0.22069017 -0.0031335442 0.027959268  1.33857649     0.18070858
## 2   1.15064473 -0.0577000699 0.486651970  1.73213460     0.08324960
## 3  -0.50633555 -0.0300891156 0.261212673 -0.93182554     0.35142669
## 4   0.21126843 -0.0551135024 0.466112384  0.39017500     0.69640715
## 5   0.12595519 -0.0095442098 0.084611240  0.46582561     0.64134034
## 6   1.30922049 -0.0488656814 0.416005238  2.10560997     0.03523825
## 7   0.70581907 -0.0507644489 0.431307151  1.15202936     0.24930901
## 8   0.13535236 -0.0007169856 0.006412863  1.69916019     0.08928900
## 9   0.47428620 -0.0165806825 0.145946648  1.28489213     0.19882998
## 10  0.97645946 -0.0337064891 0.291524842  1.87091880     0.06135634
## 11  0.68694221 -0.1236582697 0.969950667  0.82306102     0.41047329
## 12  0.93772251 -0.0449934538 0.384599458  1.58461469     0.11305391
## 13  0.06621650 -0.0026839061 0.023958142  0.44513857     0.65621962
## 14  0.64804684 -0.0341151865 0.294934838  1.25610083     0.20907942
## 15  0.06564499 -0.0005147066 0.004604571  0.97498720     0.32956661
## 16  0.25188348 -0.0370687636 0.319489333  0.51120829     0.60920521
## 17 -1.45259402 -0.0687795140 0.573277124 -1.82766069     0.06760049
## 18 -1.02774709 -0.0548725483 0.464192902 -1.42793221     0.15331137
## 19  1.00329201 -0.0588110919 0.495437677  1.50894066     0.13131395
## 20  0.24862425 -0.0189400979 0.166314737  0.65608920     0.51176674
## 21  0.28779270 -0.0332402812 0.287631352  0.59859286     0.54944442
## 22  1.14152647 -0.0687364056 0.572944337  1.59890807     0.10984103
## 23  0.68298677 -0.0116263596 0.102853220  2.16587825     0.03032049
## 24 -0.01831901 -0.0307125591 0.266453592  0.02400959     0.98084496
## 25  0.86075086 -0.0236531606 0.206702768  1.94525958     0.05174374
## 26  0.96258369 -0.0417667995 0.358224594  1.67806042     0.09333530
## 27 -1.01896533 -0.0309217943 0.268210948 -1.90782190     0.05641424
## 28  0.50571532 -0.0234450276 0.204927590  1.16892615     0.24243338
## 29  0.49265395 -0.0209596303 0.183669568  1.19844442     0.23074405
## attr(,"call")
## localmoran(x = Papua_merged$rls, listw = WLK)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1    Low-Low   Low-Low   Low-Low
## 2  High-High High-High High-High
## 3   High-Low  High-Low  High-Low
## 4   Low-High  Low-High   Low-Low
## 5   Low-High  Low-High   Low-Low
## 6    Low-Low   Low-Low   Low-Low
## 7  High-High High-High High-High
## 8    Low-Low  High-Low   Low-Low
## 9  High-High High-High High-High
## 10 High-High High-High High-High
## 11 High-High High-High High-High
## 12   Low-Low   Low-Low   Low-Low
## 13  Low-High  Low-High   Low-Low
## 14   Low-Low   Low-Low   Low-Low
## 15 High-High High-High High-High
## 16 High-High High-High High-High
## 17  High-Low  High-Low  High-Low
## 18  High-Low  High-Low  High-Low
## 19   Low-Low   Low-Low   Low-Low
## 20   Low-Low  Low-High   Low-Low
## 21  Low-High  Low-High   Low-Low
## 22   Low-Low   Low-Low   Low-Low
## 23   Low-Low   Low-Low   Low-Low
## 24 High-High High-High  High-Low
## 25 High-High High-High High-High
## 26   Low-Low   Low-Low   Low-Low
## 27  High-Low  High-Low  High-Low
## 28   Low-Low   Low-Low   Low-Low
## 29   Low-Low   Low-Low   Low-Low

Matriks bobot KNN menghasilkan p-value sebesar 0.001542

Melakukan uji Moran menggunakan matriks bobot KNN k=3

# Ekstraksi residual GWR
residual_gwr = rtg.model$SDF$gwr.e

# Membuat matriks bobot spasial berdasarkan tetangga terdekat
nb = knn2nb(knearneigh(CoordK, k=3))
listw = nb2listw(nb, style="W")

# Uji Moran's I untuk residual spasial
moran.test(residual_gwr, listw)

## 
##  Moran I test under randomisation
## 
## data:  residual_gwr  
## weights: listw    
## 
## Moran I statistic standard deviate = 1.2956, p-value = 0.09756
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12662107       -0.03571429        0.01570044

Tidak terdapat indikasi autokorealsi spasial pada residual model GWR yang berarti pola spasial telah berhasil diserap oleh model GWR

Homoskedastisitas

Uji homoskedastisitas dilakukan dengan menggunakan Uji Breusch Pagan

bptest(rtg.model$lm, weights = rtg.model$gweight)

## 
##  studentized Breusch-Pagan test
## 
## data:  rtg.model$lm
## BP = 1.4842, df = 2, p-value = 0.4761

Normalitas Residual

Pengecekan normalitas residual dilakukan dengan menggunakan uji Shapiro wilk

shapiro.test(residual_gwr)

## 
##  Shapiro-Wilk normality test
## 
## data:  residual_gwr
## W = 0.96929, p-value = 0.5406

Plot Besaran Koefisien Lokal

#Besaran koefisien lokal
koeflokal=data.frame(rtg.model$SDF)
# Menggabungkan besaran koefisien lokal dan peta Papua
Mapkoef <- cbind(Papua_sf, koeflokal)

Peta besaran koefisien lokal untuk Pegeluaran Per Kapita

ggplot(Mapkoef) +
  geom_sf(aes(fill = ppk), color = NA) +
  scale_fill_gradient(low = "yellow", high = "red",
                       limits = c(min(Mapkoef$ppk), max(Mapkoef$ppk)),
                       name = "Koefisien Lokal") +
  labs(title = "Peta Koefisien Lokal PPK") +
  theme_minimal() +
  coord_sf()

Peta besaran koefisien lokal untuk Jumlah Penduduk Miskin

ggplot(Mapkoef) +
  geom_sf(aes(fill = jpm), color = NA) +
  scale_fill_gradient(low = "yellow", high = "red",
                      limits = c(min(Mapkoef$jpm), max(Mapkoef$jpm)),
                      name = "Koefisien Lokal") +
  labs(title = "Peta Koefisien Lokal JPM ") +
  theme_minimal() +
  coord_sf()

Uji Signifikansi Lokal

t1 = rtg.model$SDF$ppk / rtg.model$SDF$ppk_se
t2 = rtg.model$SDF$jpm / rtg.model$SDF$jpm_se
id=c(1:29)

pvalue1 = 2*pt(q=abs(t1), df=24, lower.tail=FALSE) 
pvalue2 = 2*pt(q=abs(t2), df=24, lower.tail=FALSE)

Variabel x1 (Pengeluaran per Kapita)

data.frame(id, t1, pvalue1) %>% mutate(ket=ifelse(pvalue1<0.1,"Tolak H0", "Terima H0"))

##    id       t1      pvalue1       ket
## 1   1 6.934297 3.594733e-07  Tolak H0
## 2   2 4.543597 1.325581e-04  Tolak H0
## 3   3 5.296272 1.966293e-05  Tolak H0
## 4   4 7.275277 1.624151e-07  Tolak H0
## 5   5 7.109264 2.386770e-07  Tolak H0
## 6   6 8.055534 2.791670e-08  Tolak H0
## 7   7 7.553265 8.593246e-08  Tolak H0
## 8   8 8.620740 8.200044e-09  Tolak H0
## 9   9 7.189877 1.978957e-07  Tolak H0
## 10 10 7.264938 1.663386e-07  Tolak H0
## 11 11 6.678057 6.593725e-07  Tolak H0
## 12 12 8.141667 2.309891e-08  Tolak H0
## 13 13 9.341364 1.828104e-09  Tolak H0
## 14 14 8.572374 9.091250e-09  Tolak H0
## 15 15 4.361514 2.106493e-04  Tolak H0
## 16 16 1.541425 1.362965e-01 Terima H0
## 17 17 7.071095 2.608955e-07  Tolak H0
## 18 18 7.830746 4.598315e-08  Tolak H0
## 19 19 7.469229 1.040570e-07  Tolak H0
## 20 20 7.666068 6.656234e-08  Tolak H0
## 21 21 7.792970 5.003889e-08  Tolak H0
## 22 22 8.000875 3.149896e-08  Tolak H0
## 23 23 8.316830 1.576172e-08  Tolak H0
## 24 24 7.842852 4.475615e-08  Tolak H0
## 25 25 1.127176 2.708118e-01 Terima H0
## 26 26 8.814270 5.443367e-09  Tolak H0
## 27 27 8.941698 4.167413e-09  Tolak H0
## 28 28 8.129703 2.371342e-08  Tolak H0
## 29 29 8.245438 1.840952e-08  Tolak H0

Variabel x2 (Jumlah Penduduk Miskin)

data.frame(id,t2, pvalue2) %>% mutate(ket=ifelse(pvalue2<0.1,"Tolak H0", "Terima H0"))

##    id         t2      pvalue2       ket
## 1   1 -1.4287659 0.1659582169 Terima H0
## 2   2 -3.7228658 0.0010580331  Tolak H0
## 3   3 -1.1304470 0.2694589887 Terima H0
## 4   4 -1.9819902 0.0590440080  Tolak H0
## 5   5 -1.8269398 0.0801721823  Tolak H0
## 6   6 -2.7725173 0.0105824886  Tolak H0
## 7   7 -2.4957567 0.0198406411  Tolak H0
## 8   8 -2.1751799 0.0396961788  Tolak H0
## 9   9 -1.9330419 0.0651130742  Tolak H0
## 10 10 -4.6107304 0.0001117445  Tolak H0
## 11 11 -2.2739165 0.0322017358  Tolak H0
## 12 12 -2.3903388 0.0250333817  Tolak H0
## 13 13 -3.5488174 0.0016325144  Tolak H0
## 14 14 -2.4509142 0.0219139620  Tolak H0
## 15 15 -1.5781606 0.1276197970 Terima H0
## 16 16 -0.9310236 0.3611127750 Terima H0
## 17 17 -1.7243445 0.0975044731  Tolak H0
## 18 18 -2.6540004 0.0138931246  Tolak H0
## 19 19 -2.0458823 0.0518762023  Tolak H0
## 20 20 -2.3702158 0.0261569273  Tolak H0
## 21 21 -1.4584337 0.1576847518 Terima H0
## 22 22 -2.8305980 0.0092467864  Tolak H0
## 23 23 -2.9106684 0.0076654188  Tolak H0
## 24 24 -3.2576971 0.0033388884  Tolak H0
## 25 25 -0.7287260 0.4732208894 Terima H0
## 26 26 -2.7418964 0.0113581735  Tolak H0
## 27 27 -3.6701762 0.0012069544  Tolak H0
## 28 28 -1.7554987 0.0919317309  Tolak H0
## 29 29 -2.1719261 0.0399680199  Tolak H0

Peta signifikansi koefisien lokal

data1 <- data.frame(id, t1, pvalue1) %>% 
  mutate(x1 = ifelse(pvalue1 < 0.1, 1, 0))

data2 <- data.frame(id, t2, pvalue2) %>% 
  mutate(x2 = ifelse(pvalue2 < 0.1, 2, 0))

# Menggabungkan data1 dan data2 berdasarkan id
data_combined <- left_join(data1[, c("id", "x1")], data2[, c("id", "x2")], by = "id")
# Membuat kolom x3 sesuai dengan kondisi yang diberikan
data_combined <- data_combined %>%
  mutate(
    x3 = ifelse(x1 + x2 == 3, 3, 0),    # x3 bernilai 3 jika x1 + x2 = 3
    x1 = ifelse(x3 == 3, 0, x1),        # x1 diatur menjadi 0 jika x3 = 3
    x2 = ifelse(x3 == 3, 0, x2)         # x2 diatur menjadi 0 jika x3 = 3
  )

data_combined <- data_combined %>%
  mutate(Variabel = x1 + x2 + x3)

Papua_merged <- Papua_merged %>%
  left_join(data_combined, by = "id")

# Ubah variabel 'klaster' menjadi diskrit dengan 4 level
Papua_merged$variabel_diskrit <- cut(Papua_merged$Variabel, 
                                    breaks = 4, 
                                    labels = c("Tidak Signifikan", "X1", "X2", "X1 dan X2"))

# Modifikasi ggplot untuk menampilkan data diskrit
ggplot() +
  geom_sf(data = Papua_merged, aes(fill = variabel_diskrit), color = "black") +
  theme_bw() +
  scale_fill_manual(values = c("gray", "#e8f6fb", "#619fd5", "#0060ca")) +  # Warna untuk masing-masing level
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) + 
  theme(legend.position = "right",
        axis.text.x = element_blank(),  # Menghapus label sumbu x
        axis.text.y = element_blank()) +  # Menghapus label sumbu y   
  labs(title = "Signifikansi Variabel Lokal",
       fill = "Variabel yang sigfinikan")

summary(model_ols)

## 
## Call:
## lm(formula = rls ~ ppk + jpm, data = Papua_merged)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9848 -0.6989 -0.3179  1.1165  2.5734 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.1177543  1.0584450  -0.111  0.91227    
## ppk          0.0010005  0.0001136   8.809 2.78e-09 ***
## jpm         -0.0463462  0.0147302  -3.146  0.00411 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 26 degrees of freedom
## Multiple R-squared:  0.7925, Adjusted R-squared:  0.7766 
## F-statistic: 49.66 on 2 and 26 DF,  p-value: 1.32e-09

AIC(model_ols)

## [1] 111.6256

rtg.model

## Call:
## gwr(formula = rls ~ ppk + jpm, data = Papua_merged, coords = coordinates(Papua), 
##     bandwidth = bandwith_bifix, gweight = gwr.bisquare, hatmatrix = TRUE)
## Kernel function: gwr.bisquare 
## Fixed bandwidth: 3.283925 
## Summary of GWR coefficient estimates at data points:
##                     Min.     1st Qu.      Median     3rd Qu.        Max.
## X.Intercept. -2.96812856 -1.93927223 -1.08682162 -0.60104409  4.90850627
## ppk           0.00065678  0.00096750  0.00109382  0.00124482  0.00155047
## jpm          -0.13044830 -0.04527098 -0.03627370 -0.03008081 -0.02238339
##               Global
## X.Intercept. -0.1178
## ppk           0.0010
## jpm          -0.0463
## Number of data points: 29 
## Effective number of parameters (residual: 2traceS - traceS'S): 11.30889 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 17.69111 
## Sigma (residual: 2traceS - traceS'S): 1.298578 
## Effective number of parameters (model: traceS): 9.449376 
## Effective degrees of freedom (model: traceS): 19.55062 
## Sigma (model: traceS): 1.23528 
## Sigma (ML): 1.014254 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 117.6516 
## AIC (GWR p. 96, eq. 4.22): 92.56869 
## Residual sum of squares: 29.83261 
## Quasi-global R2: 0.8977016

# Mengambil prediksi dari model GWR
y_pred <- rtg.model$SDF$pred

# Mengambil nilai aktual dari variabel dependen
y_actual <- Papua_merged$rls

# Hitung Sum of Squares Total dan Sum of Squares Residual
SST <- sum((y_actual - mean(y_actual))^2)
SSR <- sum((y_actual - y_pred)^2)

# Hitung R-squared
r_squared <- 1 - (SSR / SST)
print(paste("R-squared dari model GWR:", r_squared))

## [1] "R-squared dari model GWR: 0.897701593826936"

Model GWR berhasil memenuhi seluruh asumsi, memiliki nilai AIC lebih kecil dari model OLS, dan juga memiliki nilai R-square yang lebih besar dibanding model OLS. Sehingga dapat disimpulkan bahwa model GWR berhasil memodelkan data lebih baik dibanding model OLS

Penanganan Autokorelasi Spasial pada Residual OLS

Atika Raihanah Arato

2024-10-16

Muat Package

Muat Data

Analisis Data Eksplorasi

Melihat ringkasan data secara singkat:

Membuat Peta Provinsi Papua (Sebelum Pemekaran)

Membuat Peta Sebaran RLS

Membuat Peta Sebaran PPK

Membuat Peta Sebaran JPM

Membuat Model Ordinary Least Square (OLS)

Mengecek Asumsi Model OLS

Homoskedastisitas

Multikolinearitas

Normalitas Residual

Autokorelasi Temporal

Autokorelasi Spasial

Membuat Model Geographically Weighted Regression (GWR)

Menentukan Kernel dan Bandwith yang tepat

Uji Asumsi GWR

Non Autokorelasi Spasial

Penentuan Matriks Bobot

Matriks Bobot Rook

Matriks bobot Queen

Matriks Bobot K-Nearest Neighbors

Homoskedastisitas

Normalitas Residual

Plot Besaran Koefisien Lokal

Peta besaran koefisien lokal untuk Pegeluaran Per Kapita

Peta besaran koefisien lokal untuk Jumlah Penduduk Miskin

Uji Signifikansi Lokal

Variabel x1 (Pengeluaran per Kapita)

Variabel x2 (Jumlah Penduduk Miskin)

Peta signifikansi koefisien lokal