Praktikum 11 RegSpas

Library

library(spdep)

## Warning: package 'spdep' was built under R version 4.5.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.5.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

## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE

library(sp)           
library(sf)          
library(readxl)    
library(openxlsx)

## Warning: package 'openxlsx' was built under R version 4.5.3

library(corrplot)

## corrplot 0.95 loaded

library(DescTools)

## Warning: package 'DescTools' was built under R version 4.5.3

library(nortest)     
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:DescTools':
## 
##     Recode

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.5.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

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(RColorBrewer) 
library(ggplot2) 
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ✔ readr     2.2.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ dplyr::recode() masks car::recode()
## ✖ 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(MuMIn)

## Warning: package 'MuMIn' was built under R version 4.5.3

## Registered S3 method overwritten by 'MuMIn':
##   method     from      
##   nobs.Sarlm spatialreg

library(patchwork)

Materi 1

Import Data

data <- read.csv("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/Dataset 06/data-06.csv")
head(data)

##   ID_KAB    PROPINSI              KAB     long       lat     y    x1   x2
## 1   3101 DKI JAKARTA KEPULAUAN SERIBU 106.4933 -5.796823 27.35  89.5 50.0
## 2   3171 DKI JAKARTA  JAKARTA SELATAN 106.8100 -6.272549 36.22 103.1 52.5
## 3   3172 DKI JAKARTA    JAKARTA TIMUR 106.9003 -6.255373 27.46  90.3 48.7
## 4   3173 DKI JAKARTA    JAKARTA PUSAT 106.8351 -6.181230 29.17  88.5 54.9
## 5   3174 DKI JAKARTA    JAKARTA BARAT 106.7483 -6.165469 32.16  99.6 48.1
## 6   3175 DKI JAKARTA    JAKARTA UTARA 106.8695 -6.129185 32.67 102.2 45.8

dim(data)

## [1] 119   8

peta <- st_read(dsn ="D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/Dataset 06/Jawamap", layer = "jawa")

## Reading layer `jawa' from data source 
##   `D:\Kuliah\Semester 5\Regresi Spasial\Regspas_11\Regspas_11\Dataset 06\Jawamap' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 119 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84

peta

## Simple feature collection with 119 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84
## First 10 features:
##    PROVNO KABKOTNO    PROVINSI           KABKOT ID2013
## 1      31       01 DKI JAKARTA KEPULAUAN SERIBU   3101
## 2      31       71 DKI JAKARTA  JAKARTA SELATAN   3171
## 3      31       72 DKI JAKARTA    JAKARTA TIMUR   3172
## 4      31       73 DKI JAKARTA    JAKARTA PUSAT   3173
## 5      31       74 DKI JAKARTA    JAKARTA BARAT   3174
## 6      31       75 DKI JAKARTA    JAKARTA UTARA   3175
## 7      32       01  JAWA BARAT            BOGOR   3201
## 8      32       02  JAWA BARAT         SUKABUMI   3202
## 9      32       03  JAWA BARAT          CIANJUR   3203
## 10     32       04  JAWA BARAT          BANDUNG   3204
##                          geometry
## 1  MULTIPOLYGON (((106.3842 -5...
## 2  MULTIPOLYGON (((106.8491 -6...
## 3  MULTIPOLYGON (((106.9719 -6...
## 4  MULTIPOLYGON (((106.8788 -6...
## 5  MULTIPOLYGON (((106.711 -6....
## 6  MULTIPOLYGON (((106.969 -6....
## 7  MULTIPOLYGON (((106.994 -6....
## 8  MULTIPOLYGON (((106.9652 -6...
## 9  MULTIPOLYGON (((107.2843 -6...
## 10 MULTIPOLYGON (((107.75 -6.8...

plot(peta, main = "Peta Pulau Jawa")

plot(peta[4], main = "Peta Pulau Jawa (Kabupaten)")

Eksplorasi Data

Mengecek Missing Value

data[!complete.cases(data),]

## [1] ID_KAB   PROPINSI KAB      long     lat      y        x1       x2      
## <0 rows> (or 0-length row.names)

Ringkasan Statistik Deskriptif

summary(data[,-c(1,2)])

##      KAB                 long            lat               y        
##  Length:119         Min.   :105.8   Min.   :-8.365   Min.   :13.72  
##  Class :character   1st Qu.:107.9   1st Qu.:-7.653   1st Qu.:22.07  
##  Mode  :character   Median :110.2   Median :-7.258   Median :25.40  
##                     Mean   :110.0   Mean   :-7.199   Mean   :25.53  
##                     3rd Qu.:111.8   3rd Qu.:-6.820   3rd Qu.:29.20  
##                     Max.   :114.2   Max.   :-5.797   Max.   :36.22  
##        x1               x2       
##  Min.   : 80.10   Min.   :25.10  
##  1st Qu.: 87.25   1st Qu.:33.20  
##  Median : 92.20   Median :41.20  
##  Mean   : 93.04   Mean   :40.78  
##  3rd Qu.: 99.85   3rd Qu.:49.25  
##  Max.   :105.00   Max.   :55.00

kab_max <- data %>%
  arrange(desc(y)) %>%
  slice(1)

print(paste("Kabupaten/kota dengan Y tertinggi adalah", kab_max$KAB))

## [1] "Kabupaten/kota dengan Y tertinggi adalah JAKARTA SELATAN"

kab_min <- data %>%
  arrange(y) %>%
  slice(1)

print(paste("Kabupaten/kota dengan Y terendah adalah", kab_min$KAB))

## [1] "Kabupaten/kota dengan Y terendah adalah KLATEN"

Boxplot untuk Identifikasi Pencilan

par(mfrow=c(1,3))
boxplot(data$y,main = "Sebaran Y", col="dodgerblue")
boxplot(data$x1,main = "Sebaran X1", col="dodgerblue")
boxplot(data$x2,main = "Sebaran X2", col="dodgerblue")

Hubungan Antara Peubah Respon dan Peubah Penjelas

pairs(data[,6:8], pch = 19, lower.panel = NULL)

Korelasi antar Peubah

corrplot(cor(data[,6:8]), method = "number")

cor(data[,4:7])

##             long         lat          y          x1
## long  1.00000000 -0.74541782 -0.1306153 -0.08974212
## lat  -0.74541782  1.00000000  0.1679668  0.07479095
## y    -0.13061532  0.16796681  1.0000000  0.80963262
## x1   -0.08974212  0.07479095  0.8096326  1.00000000

cor.test(data$y, data$x1)

## 
##  Pearson's product-moment correlation
## 
## data:  data$y and data$x1
## t = 14.921, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7370471 0.8637532
## sample estimates:
##       cor 
## 0.8096326

cor.test(data$y, data$x2)

## 
##  Pearson's product-moment correlation
## 
## data:  data$y and data$x2
## t = 6.5822, df = 117, p-value = 1.374e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3749256 0.6399554
## sample estimates:
##       cor 
## 0.5198398

Sebaran Spasial Data

k=101
colfunc <- colorRampPalette(c("green", "yellow","red"))

peta$y <- data$y

ggplot(peta) +
  geom_sf(aes(fill = y)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Peta Sebaran Spasial Peubah Respon Y") +
  theme_minimal()

peta$x1<- data$x1

ggplot(peta) +
  geom_sf(aes(fill = x1)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Sebaran Spasial Peubah Penjelas X1") +
  theme_minimal()

peta$x2<- data$x2

ggplot(peta) +
  geom_sf(aes(fill = x2)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Sebaran Spasial Peubah Penjelas X2") +
  theme_minimal()

Matriks Pembobot Spasial

Matriks Bobot Berdasarkan Ketetanggaan

Queen Contiguity

peta <- st_make_valid(peta)
sp.peta <- as(peta, "Spatial")

sp.polygons <- SpatialPolygons(sp.peta@polygons)
queen <- poly2nb(sp.peta, queen = T)

## Warning in poly2nb(sp.peta, queen = T): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(sp.peta, queen = T): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

queen

## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 522 
## Percentage nonzero weights: 3.68618 
## Average number of links: 4.386555 
## 1 region with no links:
## 1
## 3 disjoint connected subgraphs

W.queen <- nb2listw(queen, style='W',zero.policy=TRUE)
qc = moran.test(data$y, W.queen,randomisation=T,
alternative="greater", zero.policy = TRUE)
qc

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.queen  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.3921, p-value = 0.08195
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.085984133      -0.008547009       0.004611252

moran.mc(data$y, listw = W.queen, zero.policy = TRUE, nsim = 99)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$y 
## weights: W.queen  
## number of simulations + 1: 100 
## 
## statistic = 0.085984, observed rank = 94, p-value = 0.06
## alternative hypothesis: greater

Rook Contiguity

rook <- poly2nb(sp.peta, queen = F)

## Warning in poly2nb(sp.peta, queen = F): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(sp.peta, queen = F): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

rook

## Neighbour list object:
## Number of regions: 119 
## Number of nonzero links: 518 
## Percentage nonzero weights: 3.657934 
## Average number of links: 4.352941 
## 1 region with no links:
## 1
## 3 disjoint connected subgraphs

W.rook <- nb2listw(rook, style='W',zero.policy=TRUE)
rc = moran.test(data$y, W.rook,randomisation=T,
alternative="greater", zero.policy = TRUE)
rc

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.rook  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.5987, p-value = 0.05495
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.100243070      -0.008547009       0.004630949

moran.mc(data$y, listw = W.rook, zero.policy = TRUE, nsim = 99)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$y 
## weights: W.rook  
## number of simulations + 1: 100 
## 
## statistic = 0.10024, observed rank = 92, p-value = 0.08
## alternative hypothesis: greater

Matriks Bobot Berdasarkan Jarak

centroids <- st_centroid(peta)

## Warning: st_centroid assumes attributes are constant over geometries

longlat <- st_coordinates(centroids)

peta$long <- longlat[, 1]  
peta$lat  <- longlat[, 2]

coords <- peta[c("long","lat")]
koord <- as.data.frame(coords)
jarak<-dist(longlat, method = "euclidean")
m.jarak<-as.matrix(jarak)

K-Nearest Neighbor

W.knn<-knn2nb(knearneigh(longlat,k=5,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')

MI.knn <- moran(data$y,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(data$y,W.knn.s,randomisation = T,alternative = "greater")
mt1

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 2.1802, p-value = 0.01462
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.107628130      -0.008474576       0.002835976

Radial Distance Weight

W.rdw <-dnearneigh(longlat,0,105,longlat=TRUE)
W.rdw.s <- nb2listw(W.rdw,style='W')

MI.rdw <- moran(data$y,W.rdw.s,n=length(W.rdw.s$neighbours),S0=Szero(W.rdw.s))
mt2 = moran.test(data$y,W.rdw.s,randomisation = T,alternative = "greater")
mt2

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.rdw.s    
## 
## Moran I statistic standard deviate = 0.82869, p-value = 0.2036
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0148447999     -0.0084745763      0.0007918598

Inverse Distance Weight

alpha1=1
W.idw <-1/(m.jarak^alpha1)
diag(W.idw)<-0
rtot<-rowSums(W.idw,na.rm=TRUE)
W.idw.sd<-W.idw/rtot #row-normalized
W.idw.s = mat2listw(W.idw.sd,style='W')

MI.idw <- moran(data$y,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(data$y,W.idw.s,randomisation = T,alternative = "greater")
mt3

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 2.3687, p-value = 0.008925
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.040533196      -0.008474576       0.000428063

alpha2=2
W.idw2<-1/(m.jarak^alpha2)
diag(W.idw2)<-0
rtot<-rowSums(W.idw2,na.rm=TRUE)
W.idw2.sd<-W.idw2/rtot #row-normalized
W.idw2.s = mat2listw(W.idw2.sd,style='W')

MI.idw2 <- moran(data$y,W.idw2.s,n=length(W.idw2.s$neighbours),S0=Szero(W.idw2.s))
mt4 = moran.test(data$y,W.idw2.s,randomisation = T,alternative = "greater")
mt4

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 1.9574, p-value = 0.02515
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.088273381      -0.008474576       0.002442892

Exponential Distance Weight

alpha=1
W.edw<-exp((-alpha)*m.jarak)
diag(W.edw)<-0
rtot<-rowSums(W.edw,na.rm=TRUE)
W.edw.sd<-W.edw/rtot
W.edw.s = mat2listw(W.edw.sd,style='W')

MI.edw <- moran(data$y,W.edw.s,n=length(W.edw.s$neighbours),S0=Szero(W.edw.s))
mt5 = moran.test(data$y,W.edw.s,randomisation = T,alternative = "greater")
mt5

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 1.3022, p-value = 0.09643
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0099060556     -0.0084745763      0.0001992396

alpha2=2
W.edw2<-exp((-alpha2)*m.jarak)
diag(W.edw2)<-0
rtot<-rowSums(W.edw2,na.rm=TRUE)
W.edw2.sd<-W.edw2/rtot 
W.edw2.s = mat2listw(W.edw2.sd,style='W')
MI.edw2 <- moran(data$y,W.edw2.s,n=length(W.edw2.s$neighbours),S0=Szero(W.edw2.s))
mt6 = moran.test(data$y,W.edw2.s,randomisation = T,alternative = "greater")
mt6

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: W.edw2.s    
## 
## Moran I statistic standard deviate = 1.875, p-value = 0.03039
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0328701345     -0.0084745763      0.0004862084

Pemilihan Matriks Bobot Terbaik

MatriksBobot <- c("K-Nearest Neighbor (k=5)",
                  "Radial Distance Weight (dmax=60)",
                  "Invers Distance Weight (alpha=1)",
                  "Invers Distance Weight (alpha=2)",
                  "Exponential Distance Weight (alpha=1)",
                  "Exponential Distance Weight (alpha=2)",
                  "Rook Contiguity",
                  "Queen Contiguity")
IndeksMoran <- c(MI.knn$I,MI.rdw$I,MI.idw$I,MI.idw2$I,MI.edw$I,MI.edw2$I,rc$estimate[1],qc$estimate[1])
pv = c(mt1$p.value,mt2$p.value,mt3$p.value,mt4$p.value,mt5$p.value,mt6$p.value,rc$p.value,qc$p.value)
Matriks = cbind.data.frame(MatriksBobot, IndeksMoran,"p-value"=pv)
colnames(Matriks) <-c("Matriks Bobot", "Indeks Moran", "p-value")
Matriks

##                           Matriks Bobot Indeks Moran     p-value
## 1              K-Nearest Neighbor (k=5)  0.107628130 0.014622303
## 2      Radial Distance Weight (dmax=60)  0.014844800 0.203639555
## 3      Invers Distance Weight (alpha=1)  0.040533196 0.008925225
## 4      Invers Distance Weight (alpha=2)  0.088273381 0.025147563
## 5 Exponential Distance Weight (alpha=1)  0.009906056 0.096426630
## 6 Exponential Distance Weight (alpha=2)  0.032870135 0.030394142
## 7                       Rook Contiguity  0.100243070 0.054948820
## 8                      Queen Contiguity  0.085984133 0.081948463

Matriks[Matriks$`p-value`< 0.05,]

##                           Matriks Bobot Indeks Moran     p-value
## 1              K-Nearest Neighbor (k=5)   0.10762813 0.014622303
## 3      Invers Distance Weight (alpha=1)   0.04053320 0.008925225
## 4      Invers Distance Weight (alpha=2)   0.08827338 0.025147563
## 6 Exponential Distance Weight (alpha=2)   0.03287013 0.030394142

Model Regresi Klasik

Metode OLS

ols <- lm(y~x1+x2, data=data)
summary(ols)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5977 -0.8358  0.0369  0.9222  3.5927 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -36.66154    1.72709  -21.23   <2e-16 ***
## x1            0.54223    0.01742   31.13   <2e-16 ***
## x2            0.28804    0.01447   19.90   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.392 on 116 degrees of freedom
## Multiple R-squared:  0.922,  Adjusted R-squared:  0.9206 
## F-statistic: 685.3 on 2 and 116 DF,  p-value: < 2.2e-16

print(paste("Nilai AIC OLS", AIC(ols)))

## [1] "Nilai AIC OLS 421.320527312512"

Model Diagnostics

err.ols <- residuals(ols)
ad.test(err.ols)

## 
##  Anderson-Darling normality test
## 
## data:  err.ols
## A = 0.48312, p-value = 0.2255

hist(err.ols)

qqnorm(err.ols,datax=T)
qqline(rnorm(length(err.ols),mean(err.ols),sd(err.ols)),datax=T, col="red")

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 1.557, df = 2, p-value = 0.4591

car::vif(ols)

##      x1      x2 
## 1.00002 1.00002

Pengujian Efek Spasial

Efek Dependensi Spasial

Uji Indeks Moran

ww = W.knn.s
lm.morantest(ols, listw=ww, alternative="two.sided", zero.policy = TRUE)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = y ~ x1 + x2, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 7.7062, p-value = 1.296e-14
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.401296232     -0.008385587      0.002826244

err.ols <- residuals(ols)
merr <- moran.test(err.ols, ww,randomisation=T, 
           alternative="two.sided")
merr

## 
##  Moran I test under randomisation
## 
## data:  err.ols  
## weights: ww    
## 
## Moran I statistic standard deviate = 7.7459, p-value = 9.488e-15
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.401296232      -0.008474576       0.002798555

my <- moran.test(data$y, ww,randomisation=T, 
           alternative="two.sided")
my

## 
##  Moran I test under randomisation
## 
## data:  data$y  
## weights: ww    
## 
## Moran I statistic standard deviate = 2.1802, p-value = 0.02924
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.107628130      -0.008474576       0.002835976

moran.plot(data$y, ww, labels=data$KABKOT)

mx1 <- moran.test(data$x1, ww,randomisation=T, 
           alternative="two.sided")
mx1

## 
##  Moran I test under randomisation
## 
## data:  data$x1  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.89468, p-value = 0.371
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.039254520      -0.008474576       0.002846002

mx2 <- moran.test(data$x2, ww,randomisation=T, 
           alternative="two.sided")
mx2

## 
##  Moran I test under randomisation
## 
## data:  data$x2  
## weights: ww    
## 
## Moran I statistic standard deviate = -1.0874, p-value = 0.2768
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.066483348      -0.008474576       0.002845675

Peubah = c("Y", "X1", "X2", "Sisaan")
Indeks_Moran = c(my$estimate[1], mx1$estimate[1], mx2$estimate[1], merr$estimate[1])
p_value = c(my$p.value, mx1$p.value, mx2$p.value, merr$p.value)

df1 = cbind.data.frame(Peubah, Indeks_Moran, p_value)
colnames(df1) <- c("Peubah", "Indeks Moran", "p-value")
df1

##   Peubah Indeks Moran      p-value
## 1      Y   0.10762813 2.924461e-02
## 2     X1   0.03925452 3.709605e-01
## 3     X2  -0.06648335 2.768471e-01
## 4 Sisaan   0.40129623 9.487767e-15

Uji Lagrange Multiplier

model <- lm.LMtests(ols,listw=ww,zero.policy = TRUE, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

## Please update scripts to use lm.RStests in place of lm.LMtests

summary(model)

##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = y ~ x1 + x2, data = data)
## test weights: listw
##  
##          statistic parameter   p.value    
## RSerr      54.2452         1 1.770e-13 ***
## adjRSerr   31.9386         1 1.591e-08 ***
## RSlag      32.0571         1 1.497e-08 ***
## adjRSlag    9.7505         1  0.001793 ** 
## SARMA      63.9957         2 1.266e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Efek Heterogenitas Spasial

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 1.557, df = 2, p-value = 0.4591

Model Regresi Spasial

Model SAR

sar <- lagsarlm(y ~ x1 + x2, data = data, listw = ww)
summary(sar, Nagelkerke = T)

## 
## Call:lagsarlm(formula = y ~ x1 + x2, data = data, listw = ww)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.278572 -0.761564  0.016699  0.849727  2.695114 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -41.792688   1.851270 -22.575 < 2.2e-16
## x1            0.532067   0.015553  34.209 < 2.2e-16
## x2            0.290320   0.012806  22.671 < 2.2e-16
## 
## Rho: 0.23233, LR test value: 25.223, p-value: 5.1057e-07
## Asymptotic standard error: 0.047605
##     z-value: 4.8803, p-value: 1.0593e-06
## Wald statistic: 23.817, p-value: 1.0593e-06
## 
## Log likelihood: -194.0485 for lag model
## ML residual variance (sigma squared): 1.5135, (sigma: 1.2303)
## Nagelkerke pseudo-R-squared: 0.93687 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 398.1, (AIC for lm: 421.32)
## LM test for residual autocorrelation
## test value: 7.358, p-value: 0.0066764

err.sar<-residuals(sar)
ad.test(err.sar)

## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.28212, p-value = 0.6317

bptest.Sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 2.6875, df = 2, p-value = 0.2609

moran.test(err.sar, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: ww    
## 
## Moran I statistic standard deviate = 2.776, p-value = 0.005504
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.138466721      -0.008474576       0.002801958

Model SEM

sem<-errorsarlm(y ~ x1 + x2,data=data, listw= ww)
summary(sem)

## 
## Call:errorsarlm(formula = y ~ x1 + x2, data = data, listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.61465 -0.56031  0.16416  0.76534  2.46897 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -34.736068   1.427157 -24.339 < 2.2e-16
## x1            0.523408   0.014286  36.637 < 2.2e-16
## x2            0.283754   0.011367  24.963 < 2.2e-16
## 
## Lambda: 0.57406, LR test value: 33.057, p-value: 8.9503e-09
## Asymptotic standard error: 0.093935
##     z-value: 6.1113, p-value: 9.8842e-10
## Wald statistic: 37.348, p-value: 9.8842e-10
## 
## Log likelihood: -190.1318 for error model
## ML residual variance (sigma squared): 1.3381, (sigma: 1.1568)
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 390.26, (AIC for lm: 421.32)

err.sem<-residuals(sem)
ad.test(err.sem)

## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.75876, p-value = 0.04711

bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.89877, df = 2, p-value = 0.638

moran.test(err.sem, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.55031, p-value = 0.5821
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.037483117      -0.008474576       0.002778668

Model SARMA

SARMA <- sacsarlm(y ~ x1+x2,data=data,ww, zero.policy = TRUE)
summary(SARMA)

## 
## Call:
## sacsarlm(formula = y ~ x1 + x2, data = data, listw = ww, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.53574 -0.59853  0.10686  0.78564  2.55146 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -37.670168   2.567539 -14.672 < 2.2e-16
## x1            0.530517   0.015036  35.282 < 2.2e-16
## x2            0.286556   0.011748  24.392 < 2.2e-16
## 
## Rho: 0.083991
## Asymptotic standard error: 0.064873
##     z-value: 1.2947, p-value: 0.19543
## Lambda: 0.47951
## Asymptotic standard error: 0.11669
##     z-value: 4.1093, p-value: 3.9689e-05
## 
## LR test value: 34.119, p-value: 3.9012e-08
## 
## Log likelihood: -189.6009 for sac model
## ML residual variance (sigma squared): 1.3555, (sigma: 1.1642)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 391.2, (AIC for lm: 421.32)

err.SARMA<-residuals(SARMA)
ad.test(err.SARMA)

## 
##  Anderson-Darling normality test
## 
## data:  err.SARMA
## A = 0.70442, p-value = 0.06428

bptest.Sarlm(SARMA)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 1.2851, df = 2, p-value = 0.5259

moran.test(err.SARMA, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.36979, p-value = 0.7115
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.027983670      -0.008474576       0.002783258

Model SDM

sdm <- lagsarlm(y ~ x1 + x2, data=data, ww, zero.policy = TRUE, type="mixed");
summary(sdm)

## 
## Call:lagsarlm(formula = y ~ x1 + x2, data = data, listw = ww, type = "mixed", 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -4.460816 -0.587583  0.063713  0.803323  2.489449 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -23.011394   5.406184 -4.2565 2.077e-05
## x1            0.528068   0.014497 36.4268 < 2.2e-16
## x2            0.287681   0.012170 23.6394 < 2.2e-16
## lag.x1       -0.226018   0.064788 -3.4886 0.0004856
## lag.x2       -0.113673   0.045416 -2.5029 0.0123175
## 
## Rho: 0.52166, LR test value: 25.205, p-value: 5.1548e-07
## Asymptotic standard error: 0.10037
##     z-value: 5.1973, p-value: 2.0223e-07
## Wald statistic: 27.012, p-value: 2.0223e-07
## 
## Log likelihood: -188.6373 for mixed model
## ML residual variance (sigma squared): 1.3227, (sigma: 1.1501)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 391.27, (AIC for lm: 414.48)
## LM test for residual autocorrelation
## test value: 5.7604, p-value: 0.016392

err.sdm<-sdm$residuals 
ad.test(err.sdm)

## 
##  Anderson-Darling normality test
## 
## data:  err.sdm
## A = 0.64985, p-value = 0.08783

hist(err.sdm)

qqnorm(err.sdm,datax=T) 
qqline(rnorm(length(err.sdm),mean(err.sdm),sd(err.sdm)),datax=T, col="blue")

bptest.Sarlm(sdm)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 2.774, df = 4, p-value = 0.5963

moran.test(err.sdm, listw=ww, alternative="two.sided", zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  err.sdm  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.59546, p-value = 0.5515
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.039895217      -0.008474576       0.002784333

Goodness of Fit (Model Terbaik)

df <- data.frame("Model" = c("OLS (Regresi KlasiK)","SAR","SEM","SARMA", "SDM"),
           "AIC" = c(AIC(ols),AIC(sar), AIC(sem),AIC(SARMA), AIC(sdm)))

df

##                  Model      AIC
## 1 OLS (Regresi KlasiK) 421.3205
## 2                  SAR 398.0970
## 3                  SEM 390.2637
## 4                SARMA 391.2017
## 5                  SDM 391.2745

Interpretasi dan Kesimpulan

sum <- summary(SARMA)
sum

## 
## Call:
## sacsarlm(formula = y ~ x1 + x2, data = data, listw = ww, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.53574 -0.59853  0.10686  0.78564  2.55146 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -37.670168   2.567539 -14.672 < 2.2e-16
## x1            0.530517   0.015036  35.282 < 2.2e-16
## x2            0.286556   0.011748  24.392 < 2.2e-16
## 
## Rho: 0.083991
## Asymptotic standard error: 0.064873
##     z-value: 1.2947, p-value: 0.19543
## Lambda: 0.47951
## Asymptotic standard error: 0.11669
##     z-value: 4.1093, p-value: 3.9689e-05
## 
## LR test value: 34.119, p-value: 3.9012e-08
## 
## Log likelihood: -189.6009 for sac model
## ML residual variance (sigma squared): 1.3555, (sigma: 1.1642)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 391.2, (AIC for lm: 421.32)

Im = impacts(SARMA, listw = ww)
Im

## Impact measures (sac, exact):
##             Direct   Indirect     Total
## x1 dy/dx 0.5311153 0.04804611 0.5791614
## x2 dy/dx 0.2868790 0.02595184 0.3128308

koef = sum$coefficients[-1]
diref = Im$direct
umbal = diref-koef
cbind.data.frame(Koefisien=koef,EfekLangsung=diref,UmpanBalik=umbal)

##    Koefisien EfekLangsung   UmpanBalik
## x1 0.5305173    0.5311153 0.0005980484
## x2 0.2865560    0.2868790 0.0003230325

Peubah Penjelas X1

Pengaruh langsung dari peubah X1 adalah sebesar 0.531, artinya Jika rata-rata X1 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah tersebut akan meningkat sebesar 0.531, jika peubah lainnya tetap.

Pengaruh tidak langsung dari peubah X1 bernilai 0.049, artinya jika Jika rata-rata X1 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah-j (wilayah yang berbeda) akan meningkat sebesar 0.049, jika peubah lainnya tetap.

Pengaruh total dari peubah X1 bernilai 0.580, artinya apabila peningkatan X1 terjadi di seluruh wilayah maka akan meningkatkan Y diseluruh wilayah dengan rata-rata peningkatan sebesar 0.580. Efek umpan baliknya sebesar 0.0006.

Peubah Penjelas X2

Pengaruh langsung dari peubah X2 adalah sebesar 0.287, artinya Jika rata-rata X2 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah tersebut akan meningkat sebesar 0.287, jika peubah lainnya tetap.

Pengaruh tidak langsung dari X2 tersebut bernilai 0.026, artinya jika Jika rata-rata X2 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah-j (wilayah yang berbeda) akan meningkat sebesar 0.026, jika peubah lainnya tetap.

Pengaruh total dari peubah X2 adalah sebesar 0.313, artinya apabila peningkatan X2 terjadi di seluruh wilayah maka akan meningkatkan Y diseluruh wilayah dengan rata-rata peningkatan sebesar 0.313, jika peubah lainnya tetap. Efek umpan baliknya sebesar 0.000341.

Im2 = impacts(sar, listw = ww)
Im2

## Impact measures (lag, exact):
##             Direct   Indirect     Total
## x1 dy/dx 0.5370876 0.15600347 0.6930911
## x2 dy/dx 0.2930588 0.08512239 0.3781812

Im3 = impacts(sdm, listw = ww)
Im3

## Impact measures (mixed, exact):
##             Direct   Indirect     Total
## x1 dy/dx 0.5339407 0.09751864 0.6314593
## x2 dy/dx 0.2920035 0.07177461 0.3637781

Materi 2

Model Prediksi GRP Data China

Persiapan Data

shp_china <- st_read("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/data china/China29/China29.shp")

## Reading layer `China29' from data source 
##   `D:\Kuliah\Semester 5\Regresi Spasial\Regspas_11\Regspas_11\data china\China29\China29.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 29 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 73.5577 ymin: 20.22264 xmax: 134.7739 ymax: 53.56086
## Geodetic CRS:  WGS 84

data_china <- read_excel("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/data china/Data_China.xlsx")

china_join <- shp_china %>%
  left_join(data_china,
            by = c("NAME_1" = "Region"))
head(china_join)

## Simple feature collection with 6 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 92.76416 ymin: 20.22264 xmax: 120.8376 ymax: 42.77393
## Geodetic CRS:  WGS 84
##   ID_0 ISO NAME_0 ID_1    NAME_1 HASC_1 CCN_1 CCA_1    TYPE_1    ENGTYPE_1
## 1   49 CHN  China    1     Anhui  CN.AH     0  <NA>     Sh?ng     Province
## 2   49 CHN  China    2   Beijing  CN.BJ     0  <NA> Zhíxiáshì Municipality
## 3   49 CHN  China    3 Chongqing  CN.CQ     0  <NA> Zhíxiáshì Municipality
## 4   49 CHN  China    4    Fujian  CN.FJ     0  <NA>     Sh?ng     Province
## 5   49 CHN  China    5     Gansu  CN.GS     0  <NA>     Sh?ng     Province
## 6   49 CHN  China    6 Guangdong  CN.GD     0  <NA>     Sh?ng     Province
##   NL_NAME_1 VARNAME_1     Long      Lat   GRP Level_of_Urbanization
## 1     ??|??     ?nhu? 117.2262 31.82579 28792                 46.50
## 2     ??|??   B?ij?ng 116.4107 40.18491 87475                 86.20
## 3     ??|?? Chóngqìng 107.8748 30.05865 38914                 56.98
## 4        ??    Fújiàn 117.9815 26.08641 52763                 59.60
## 5     ??|??     G?nsù 100.9339 37.82079 21978                 38.75
## 6     ??|?? Gu?ngd?ng 113.4164 23.34225 54095                 67.40
##                         geometry
## 1 MULTIPOLYGON (((116.4263 34...
## 2 MULTIPOLYGON (((116.6669 40...
## 3 MULTIPOLYGON (((108.5419 32...
## 4 MULTIPOLYGON (((117.689 23....
## 5 MULTIPOLYGON (((97.18472 42...
## 6 MULTIPOLYGON (((110.1182 20...

Eksplorasi Data

GRP <- ggplot(data = china_join) +
  geom_sf(
    aes(fill = GRP),
    color = "white",
    linewidth = 0.2
  ) +
  
  scale_fill_viridis_c(
    option = "magma",
    name = "GRP"
  ) +
  
  labs(
    title = "Gross Regional Product",
    subtitle = "Data kontinu per provinsi",
    x = "Longitude",
    y = "Latitude"
  ) +
  
  theme_minimal() +
  
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 11)
  )

Urbanization <- ggplot(data = china_join) +
  geom_sf(
    aes(fill = `Level_of_Urbanization`),
    color = "white",
    linewidth = 0.2
  ) +
  
  scale_fill_viridis_c(
    option = "magma",
    name = "Level of Urbanization"
  ) +
  
  labs(
    title = "Level of Urbanization",
    subtitle = "Data kontinu per provinsi",
    x = "Longitude",
    y = "Latitude"
  ) +
  
  theme_minimal() +
  
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 11)
  )

GRP | Urbanization

str(china_join)

## Classes 'sf' and 'data.frame':   29 obs. of  17 variables:
##  $ ID_0                 : chr  "49" "49" "49" "49" ...
##  $ ISO                  : chr  "CHN" "CHN" "CHN" "CHN" ...
##  $ NAME_0               : chr  "China" "China" "China" "China" ...
##  $ ID_1                 : chr  "1" "2" "3" "4" ...
##  $ NAME_1               : chr  "Anhui" "Beijing" "Chongqing" "Fujian" ...
##  $ HASC_1               : chr  "CN.AH" "CN.BJ" "CN.CQ" "CN.FJ" ...
##  $ CCN_1                : chr  "0" "0" "0" "0" ...
##  $ CCA_1                : chr  NA NA NA NA ...
##  $ TYPE_1               : chr  "Sh?ng" "Zhíxiáshì" "Zhíxiáshì" "Sh?ng" ...
##  $ ENGTYPE_1            : chr  "Province" "Municipality" "Municipality" "Province" ...
##  $ NL_NAME_1            : chr  "??|??" "??|??" "??|??" "??" ...
##  $ VARNAME_1            : chr  "?nhu?" "B?ij?ng" "Chóngqìng" "Fújiàn" ...
##  $ Long                 : num  117 116 108 118 101 ...
##  $ Lat                  : num  31.8 40.2 30.1 26.1 37.8 ...
##  $ GRP                  : num  28792 87475 38914 52763 21978 ...
##  $ Level_of_Urbanization: num  46.5 86.2 57 59.6 38.8 ...
##  $ geometry             :sfc_MULTIPOLYGON of length 29; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:7291, 1:2] 116 116 116 116 116 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:16] "ID_0" "ISO" "NAME_0" "ID_1" ...

colSums(is.na(china_join))

##                  ID_0                   ISO                NAME_0 
##                     0                     0                     0 
##                  ID_1                NAME_1                HASC_1 
##                     0                     0                     0 
##                 CCN_1                 CCA_1                TYPE_1 
##                     0                    29                     0 
##             ENGTYPE_1             NL_NAME_1             VARNAME_1 
##                     0                     0                     0 
##                  Long                   Lat                   GRP 
##                     0                     0                     0 
## Level_of_Urbanization              geometry 
##                     0                     0

Regresi OLS

china.lm<- lm(GRP ~ Level_of_Urbanization,data=china_join)
summary(china.lm)

## 
## Call:
## lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12050.9  -4410.4   -374.8   2961.6  14948.1 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -31895.5     5197.6  -6.137 1.48e-06 ***
## Level_of_Urbanization   1400.0       92.6  15.118 1.06e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6595 on 27 degrees of freedom
## Multiple R-squared:  0.8944, Adjusted R-squared:  0.8904 
## F-statistic: 228.6 on 1 and 27 DF,  p-value: 1.064e-14

Matriks Pembobot Spasial

rail <- read_excel("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/data china/Data_China.xlsx", sheet = "Railway Distance", skip = 1)

## New names:
## • `` -> `...1`

rail_matrix <- as.matrix(rail[, 3:31])
rownames(rail_matrix) <- colnames(rail_matrix)

W_rail <- 1 / rail_matrix
diag(W_rail) <- 0

W_rail_std <- W_rail / rowSums(W_rail)

W <- mat2listw(W_rail_std, style = "W")

Uji Morans’I

col.moran <- lm.morantest(china.lm, W)
col.moran

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
## weights: W
## 
## Moran I statistic standard deviate = 1.1234, p-value = 0.1306
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.010629595     -0.033552435      0.001546711

Uji Lagrange Multiplier

china.lagrange1 <- lm.LMtests(china.lm, W, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

## Please update scripts to use lm.RStests in place of lm.LMtests

summary(china.lagrange1)

##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = GRP ~ Level_of_Urbanization, data = china_join)
## test weights: listw
##  
##          statistic parameter p.value
## RSerr    0.0277967         1  0.8676
## adjRSerr 0.0726855         1  0.7875
## RSlag    0.0035418         1  0.9525
## adjRSlag 0.0484307         1  0.8258
## SARMA    0.0762273         2  0.9626

SLX Model

china.SLX <- lmSLX(GRP ~ Level_of_Urbanization, data = china_join, listw =  W, Durbin = TRUE);
summary(china.SLX)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                            Estimate    Std. Error  t value     Pr(>|t|)  
## (Intercept)                -2.335e+04   4.130e+04  -5.654e-01   5.766e-01
## Level_of_Urbanization       1.387e+03   1.139e+02   1.217e+01   3.074e-12
## lag.Level_of_Urbanization  -1.418e+02   6.801e+02  -2.086e-01   8.364e-01

LR Test: SLX vs OLS Regression

logLik(china.SLX)

## 'log Lik.' -295.1173 (df=4)

logLik(china.lm)

## 'log Lik.' -295.1415 (df=3)

LR <- 2 * (as.numeric(logLik(china.SLX)) - 
           as.numeric(logLik(china.lm)))

df <- attr(logLik(china.SLX), "df") - 
      attr(logLik(china.lm), "df")

pchisq(LR, df = df, lower.tail = FALSE)

## [1] 0.8257453

AICc(china.lm)

## [1] 597.2431

AICc(china.SLX)

## [1] 599.9012

Model Dependensi Spasial (Jawa Barat)

Data Preprocessing

shp_jabar <- st_read("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/JABAR/Jabar2/Jabar2.shp")

## Reading layer `Jabar2' from data source 
##   `D:\Kuliah\Semester 5\Regresi Spasial\Regspas_11\Regspas_11\JABAR\Jabar2\Jabar2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS:           NA

data_jabar <- read_excel("D:/Kuliah/Semester 5/Regresi Spasial/Regspas_11/Regspas_11/JABAR/data_jabar.xlsx")

data_jabar$KABKOTNO <- as.character(data_jabar$KABKOTNO)
jabar <- shp_jabar %>%
  left_join(data_jabar, by = "KABKOTNO") %>%
  na.omit()

colnames(jabar)

##  [1] "PROVNO.x"   "KABKOTNO"   "KODE2010.x" "PROVINSI.x" "KABKOT.x"  
##  [6] "SUMBER"     "IDSP2010.x" "PROVNO.y"   "KODE2010.y" "PROVINSI.y"
## [11] "KABKOT.y"   "IDSP2010.y" "Long"       "Lat"        "p.miskin15"
## [16] "p.miskin16" "j.miskin15" "j.miskin16" "AHH2015"    "AHH2016"   
## [21] "EYS2015"    "EYS2016"    "MYS2015"    "MYS2016"    "EXP2015"   
## [26] "EXP2016"    "APM.SD15"   "APM.SMP15"  "APM.SMA15"  "APM.PT15"  
## [31] "APK.SD15"   "APK.SMP15"  "APK.SMA15"  "APK.PT15"   "APS.USIA15"
## [36] "APS.USIA2"  "APS.USIA3"  "APS.USIA4"  "geometry"

Exploratory Analysis

plot(jabar[,c(
  "p.miskin15",
  "AHH2015",
  "EYS2015",
  "MYS2015",
  "EXP2015",
  "APM.SD15",
  "APM.SMP15",
  "APM.SMA15",
  "APM.PT15",
  "APK.SD15",
  "APK.SMP15",
  "APK.SMA15",
  "APK.PT15",
  "APS.USIA15",
  "APS.USIA2",
  "APS.USIA3",
  "APS.USIA4"
)])

## Warning: plotting the first 9 out of 17 attributes; use max.plot = 17 to plot
## all

Matriks Pembobot Spasial

nb_jabar <- suppressWarnings(poly2nb(jabar))

lw_jabar <- nb2listw(
  nb_jabar,
  style = "W",
  zero.policy = TRUE
)

Penentuan Peubah (Indeks Moran)

variabel_list <- c(
  "p.miskin15",
  "AHH2015",
  "EYS2015",
  "MYS2015",
  "EXP2015",
  "APM.SD15",
  "APM.SMP15",
  "APM.SMA15",
  "APM.PT15",
  "APK.SD15",
  "APK.SMP15",
  "APK.SMA15",
  "APK.PT15",
  "APS.USIA15",
  "APS.USIA2",
  "APS.USIA3",
  "APS.USIA4"
)

variabel_Y <- "p.miskin15"

variabel_X <- setdiff(
  variabel_list,
  variabel_Y
)

hasil_model <- data.frame(
  Formula = character(),
  Jumlah_X = numeric(),
  Moran_I = numeric(),
  P_Value = numeric(),
  stringsAsFactors = FALSE
)

maksimum_X <- 6

for (k in 1:maksimum_X) {
  kombinasi_x <- combn(
    variabel_X,
    k,
    simplify = FALSE
  )
  
  for (x_set in kombinasi_x) {
    
    tryCatch({
      formula_text <- paste(
        variabel_Y,
        "~",
        paste(x_set, collapse = " + ")
      )
      
      rumus <- as.formula(formula_text)
      
      model_ols <- lm(
        rumus,
        data = jabar
      )
      
      uji_moran <- lm.morantest(
        model_ols,
        lw_jabar,
        zero.policy = TRUE
      )
      
      moran_i <- as.numeric(
        uji_moran$estimate["Observed Moran I"]
      )
      
      p_val <- uji_moran$p.value
      
      if (p_val < 0.05) {
        
        hasil_model <- rbind(
          hasil_model,
          
          data.frame(
            Formula = formula_text,
            Jumlah_X = k,
            Moran_I = moran_i,
            P_Value = p_val
          )
        )
      }
      
    }, error = function(e) {
      
      NULL
      
    })
  }
}

hasil_model <- hasil_model %>%
  arrange(P_Value)


if (nrow(hasil_model) == 0) {
  
  cat("Tidak ada model signifikan.\n")
  
} else {
  
  model_terbaik <- hasil_model %>%
    slice(1)
  
  cat(
    "Formula              : ",
    as.character(model_terbaik$Formula),
    "\n",
    sep = ""
  )
  
  cat(
    "Jumlah Variabel X    : ",
    model_terbaik$Jumlah_X,
    "\n",
    sep = ""
  )
  
  cat(
    "Moran's I            : ",
    round(model_terbaik$Moran_I, 5),
    "\n",
    sep = ""
  )
  
  cat(
    "P-Value              : ",
    format(
      model_terbaik$P_Value,
      scientific = FALSE
    ),
    "\n",
    sep = ""
  )
}

## Formula              : p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15
## Jumlah Variabel X    : 5
## Moran's I            : 0.5012
## P-Value              : 0.005018003

OLS Regression

jabar.lm<-lm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15, data=jabar)
summary(jabar.lm)

## 
## Call:
## lm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + 
##     APS.USIA15, data = jabar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2541 -1.5129  0.7468  1.7253  2.9766 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -263.5186   236.1804  -1.116   0.2883  
## EYS2015        2.7220     2.1113   1.289   0.2238  
## APM.SD15      -0.5289     0.6939  -0.762   0.4619  
## APK.SD15       0.4526     0.2946   1.536   0.1527  
## APK.PT15      -0.3228     0.1731  -1.865   0.0891 .
## APS.USIA15     2.4849     2.5164   0.987   0.3446  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.114 on 11 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.3816 
## F-statistic: 2.975 on 5 and 11 DF,  p-value: 0.06127

Uji Lagrange Multiplier

jabar.lagrange <- lm.LMtests(jabar.lm, lw_jabar, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

## Please update scripts to use lm.RStests in place of lm.LMtests

summary(jabar.lagrange)

##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 +
## APK.PT15 + APS.USIA15, data = jabar)
## test weights: listw
##  
##          statistic parameter  p.value   
## RSerr      7.23955         1 0.007132 **
## adjRSerr   5.48939         1 0.019132 * 
## RSlag      2.01129         1 0.156133   
## adjRSlag   0.26113         1 0.609346   
## SARMA      7.50067         2 0.023510 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Spatial Error Model (SEM)

jabar.err <- errorsarlm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar,data=jabar, zero.policy = TRUE)
summary(jabar.err)

## 
## Call:errorsarlm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + 
##     APK.PT15 + APS.USIA15, data = jabar, listw = lw_jabar, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.96720 -0.97706  0.17825  0.76731  1.45103 
## 
## Type: error 
## Regions with no neighbours included:
##  18 19 21 25 26 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -425.458866   84.564213 -5.0312 4.874e-07
## EYS2015        4.311729    1.035451  4.1641 3.126e-05
## APM.SD15      -1.521207    0.226462 -6.7173 1.852e-11
## APK.SD15       0.640932    0.101930  6.2880 3.216e-10
## APK.PT15      -0.358322    0.069408 -5.1626 2.436e-07
## APS.USIA15     4.680477    0.828208  5.6513 1.592e-08
## 
## Lambda: 0.9059, LR test value: 23.088, p-value: 1.5476e-06
## Asymptotic standard error: 0.060606
##     z-value: 14.947, p-value: < 2.22e-16
## Wald statistic: 223.43, p-value: < 2.22e-16
## 
## Log likelihood: -28.18863 for error model
## ML residual variance (sigma squared): 1.0473, (sigma: 1.0234)
## Number of observations: 17 
## Number of parameters estimated: 8 
## AIC: 72.377, (AIC for lm: 93.465)

Model SLX

jabar.SLX <- lmSLX(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar,data=jabar, zero.policy = TRUE, Durbin = TRUE);
summary(jabar.SLX)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                 Estimate    Std. Error  t value     Pr(>|t|)  
## (Intercept)     -145.59982   307.35728    -0.47372     0.65245
## EYS2015            2.63817     4.59362     0.57431     0.58663
## APM.SD15           0.26488     0.97666     0.27122     0.79532
## APK.SD15           0.59059     0.23882     2.47293     0.04827
## APK.PT15          -0.20577     0.38277    -0.53759     0.61021
## APS.USIA15         0.36388     3.17195     0.11472     0.91241
## lag.EYS2015       -4.02200     9.31944    -0.43157     0.68113
## lag.APM.SD15       3.42980     1.21392     2.82538     0.03014
## lag.APK.SD15      -0.41090     0.69159    -0.59414     0.57412
## lag.APK.PT15       0.38003     0.82137     0.46268     0.65990
## lag.APS.USIA15    -2.48727     1.03245    -2.40910     0.05264

Model SARMA

jabar.sarma <- sacsarlm(p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + APK.PT15 + APS.USIA15,lw_jabar, data=jabar, zero.policy = TRUE)
summary(jabar.sarma)

## 
## Call:sacsarlm(formula = p.miskin15 ~ EYS2015 + APM.SD15 + APK.SD15 + 
##     APK.PT15 + APS.USIA15, data = jabar, listw = lw_jabar, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.837342 -0.900247 -0.050073  0.685523  1.606525 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -442.074286   83.148799 -5.3167 1.057e-07
## EYS2015        4.793132    1.140231  4.2037 2.626e-05
## APM.SD15      -1.707710    0.293038 -5.8276 5.623e-09
## APK.SD15       0.672837    0.103147  6.5231 6.886e-11
## APK.PT15      -0.392055    0.076727 -5.1097 3.226e-07
## APS.USIA15     4.939024    0.838588  5.8897 3.869e-09
## 
## Rho: 0.20377
## Asymptotic standard error: 0.19972
##     z-value: 1.0202, p-value: 0.30761
## Lambda: 0.91256
## Asymptotic standard error: 0.065071
##     z-value: 14.024, p-value: < 2.22e-16
## 
## LR test value: 23.918, p-value: 6.4015e-06
## 
## Log likelihood: -27.77358 for sac model
## ML residual variance (sigma squared): 0.97176, (sigma: 0.98578)
## Number of observations: 17 
## Number of parameters estimated: 9 
## AIC: 73.547, (AIC for lm: 93.465)

LR Test: SAR vs SARMA

anova(jabar.err,jabar.sarma)

##             Model df    AIC  logLik Test L.Ratio p-value
## jabar.err       1  8 72.377 -28.189    1                
## jabar.sarma     2  9 73.547 -27.774    2 0.83008 0.36225

Akaike Information Criterion Corrected (AICc)

AICc(jabar.lm)

## [1] 105.9096

AICc(jabar.err)

## [1] 90.37725

AICc(jabar.sarma)

## [1] 99.26145