Regspas kalimatnan

Set up library

library(spdep)        # Analisis data spasial: spatial autocorrelation, spatial regression

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

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.3.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(sp)           # Struktur data untuk manipulasi data spasial (poligon, titik, dll)

## Warning: package 'sp' was built under R version 4.3.3

library(sf)           # Mengelola data spasial dengan standar Simple Features (sf)
library(readxl)       # Membaca file Excel (format .xls dan .xlsx)

## Warning: package 'readxl' was built under R version 4.3.2

library(openxlsx)     # Membuat, membaca, dan menulis file Excel (lebih cepat dari readxl)

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

library(corrplot)     # Membuat visualisasi korelasi dalam bentuk plot

## Warning: package 'corrplot' was built under R version 4.3.2

## corrplot 0.92 loaded

library(DescTools)    # Kumpulan fungsi statistik deskriptif dan inferensial

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

library(nortest)      # Melakukan uji normalitas (Anderson-Darling, Lilliefors, dll)
library(car)          # Fungsi regresi linear dan diagnostik model (vif, durbinWatsonTest, dll)

## Warning: package 'car' was built under R version 4.3.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.2

## 
## Attaching package: 'car'

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

library(spatialreg)   # Estimasi model regresi spasial (SAR, SEM, SLM)

## Warning: package 'spatialreg' was built under R version 4.3.3

## Loading required package: Matrix

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(lmtest)       # Fungsi untuk menguji model regresi (heteroskedastisitas, serial correlation)

## Warning: package 'lmtest' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.2

## 
## Attaching package: 'zoo'

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

library(dplyr)        # Manipulasi data seperti filter, select, mutate, dan summarise

## Warning: package 'dplyr' was built under R version 4.3.2

## 
## 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) # Skema palet warna yang dapat digunakan dalam visualisasi
library(ggplot2)      # Membuat grafik statistik dengan pendekatan Grammar of Graphics

## Warning: package 'ggplot2' was built under R version 4.3.3

Import data

Data Frame

data <- read_excel("C:/Users/Admin/Downloads/Data Regspas Kalimantan UAS.xlsx")
data

## # A tibble: 56 × 8
##    Daerah        IPM    X1    X2    X3    X4    X5    X6
##    <chr>       <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 Sambas       68.7  7.08  6.75 11057  74.0  12.7  5.04
##  2 Bengkayang   69.5  6.28  7.22  9859  74.2  12.2  2.92
##  3 Landak       68.1  9.97  7.24  8423  74.0  12.5  2.24
##  4 Mempawah     67.9  5.21  7.2   8690  74.0  12.9  7.33
##  5 Sanggau      67.8  4.79  7.44  9213  74.0  11.9  3.86
##  6 Ketapang     68.7  9.25  7.55  9984  73.6  12.0  6.57
##  7 Sintang      68.7  8.18  7.64  9128  74.1  12.3  2.92
##  8 Kapuas Hulu  67.7  8.16  7.82  8055  73.3  12.2  2.19
##  9 Sekadau      66.3  5.9   7.13  7981  73.9  11.9  2.29
## 10 Melawi       67.8 11.1   7.42  9155  73.7  11.3  2.46
## # ℹ 46 more rows

SHP

shp<- ("C:/Users/Admin/Downloads/Kalimantan Tengah/SHP_Kalimantan.shp")
shp<- st_read(shp)

## Reading layer `SHP_Kalimantan' from data source 
##   `C:\Users\Admin\Downloads\Kalimantan Tengah\SHP_Kalimantan.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 56 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6778 ymin: -4.744813 xmax: 118.9891 ymax: 4.408303
## Geodetic CRS:  WGS 84

head(shp)

## Simple feature collection with 6 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 114.3466 ymin: -3.730723 xmax: 115.8575 ymax: -0.0315456
## Geodetic CRS:  WGS 84
##   Shape_Leng Shape_Area        ADM2_EN ADM2_PCODE ADM2_REF ADM2ALT1EN
## 1   2.106946  0.1546742       Balangan     ID6311     <NA>       <NA>
## 2   4.691104  0.3845014         Banjar     ID6303     <NA>       <NA>
## 3   3.058313  0.1936082   Barito Kuala     ID6304     <NA>       <NA>
## 4   4.518498  0.5392015 Barito Selatan     ID6204     <NA>       <NA>
## 5   2.532642  0.2676616   Barito Timur     ID6212     <NA>       <NA>
## 6   5.214343  0.8071958   Barito Utara     ID6205     <NA>       <NA>
##   ADM2ALT2EN            ADM1_EN ADM1_PCODE   ADM0_EN ADM0_PCODE       date
## 1       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 2       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 3       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 4       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
## 5       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
## 6       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
##      validOn    validTo                       geometry
## 1 2020-04-01 -001-11-30 MULTIPOLYGON (((115.7292 -1...
## 2 2020-04-01 -001-11-30 MULTIPOLYGON (((115.5562 -2...
## 3 2020-04-01 -001-11-30 MULTIPOLYGON (((114.4129 -3...
## 4 2020-04-01 -001-11-30 MULTIPOLYGON (((115.0968 -1...
## 5 2020-04-01 -001-11-30 MULTIPOLYGON (((115.4157 -1...
## 6 2020-04-01 -001-11-30 MULTIPOLYGON (((115.3015 -0...

Merged Data

head(shp)

## Simple feature collection with 6 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 114.3466 ymin: -3.730723 xmax: 115.8575 ymax: -0.0315456
## Geodetic CRS:  WGS 84
##   Shape_Leng Shape_Area        ADM2_EN ADM2_PCODE ADM2_REF ADM2ALT1EN
## 1   2.106946  0.1546742       Balangan     ID6311     <NA>       <NA>
## 2   4.691104  0.3845014         Banjar     ID6303     <NA>       <NA>
## 3   3.058313  0.1936082   Barito Kuala     ID6304     <NA>       <NA>
## 4   4.518498  0.5392015 Barito Selatan     ID6204     <NA>       <NA>
## 5   2.532642  0.2676616   Barito Timur     ID6212     <NA>       <NA>
## 6   5.214343  0.8071958   Barito Utara     ID6205     <NA>       <NA>
##   ADM2ALT2EN            ADM1_EN ADM1_PCODE   ADM0_EN ADM0_PCODE       date
## 1       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 2       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 3       <NA> Kalimantan Selatan       ID63 Indonesia         ID 2019-12-20
## 4       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
## 5       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
## 6       <NA>  Kalimantan Tengah       ID62 Indonesia         ID 2019-12-20
##      validOn    validTo                       geometry
## 1 2020-04-01 -001-11-30 MULTIPOLYGON (((115.7292 -1...
## 2 2020-04-01 -001-11-30 MULTIPOLYGON (((115.5562 -2...
## 3 2020-04-01 -001-11-30 MULTIPOLYGON (((114.4129 -3...
## 4 2020-04-01 -001-11-30 MULTIPOLYGON (((115.0968 -1...
## 5 2020-04-01 -001-11-30 MULTIPOLYGON (((115.4157 -1...
## 6 2020-04-01 -001-11-30 MULTIPOLYGON (((115.3015 -0...

mapauas <- full_join(data, shp, by = c("Daerah" = "ADM2_EN"))
mapauas

## # A tibble: 56 × 22
##    Daerah        IPM    X1    X2    X3    X4    X5    X6 Shape_Leng Shape_Area
##    <chr>       <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>      <dbl>      <dbl>
##  1 Sambas       68.7  7.08  6.75 11057  74.0  12.7  5.04       4.06      0.473
##  2 Bengkayang   69.5  6.28  7.22  9859  74.2  12.2  2.92       5.40      0.458
##  3 Landak       68.1  9.97  7.24  8423  74.0  12.5  2.24       5.27      0.670
##  4 Mempawah     67.9  5.21  7.2   8690  74.0  12.9  7.33       2.66      0.163
##  5 Sanggau      67.8  4.79  7.44  9213  74.0  11.9  3.86       7.10      1.04 
##  6 Ketapang     68.7  9.25  7.55  9984  73.6  12.0  6.57      12.0       2.42 
##  7 Sintang      68.7  8.18  7.64  9128  74.1  12.3  2.92      10.8       1.82 
##  8 Kapuas Hulu  67.7  8.16  7.82  8055  73.3  12.2  2.19       9.78      2.53 
##  9 Sekadau      66.3  5.9   7.13  7981  73.9  11.9  2.29       4.17      0.463
## 10 Melawi       67.8 11.1   7.42  9155  73.7  11.3  2.46       5.12      0.814
## # ℹ 46 more rows
## # ℹ 12 more variables: ADM2_PCODE <chr>, ADM2_REF <chr>, ADM2ALT1EN <chr>,
## #   ADM2ALT2EN <chr>, ADM1_EN <chr>, ADM1_PCODE <chr>, ADM0_EN <chr>,
## #   ADM0_PCODE <chr>, date <date>, validOn <date>, validTo <date>,
## #   geometry <MULTIPOLYGON [°]>

Eksplorasi Data

Box Plot

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

boxplot(data$X3,main = "Sebaran X1", col="dodgerblue")
boxplot(data$X4,main = "Sebaran X2", col="dodgerblue")
boxplot(data$X5,main = "Sebaran X1", col="dodgerblue")

boxplot(data$X6,main = "Sebaran X2", col="dodgerblue")

IPM-nya ada pencilanya trus jaraknya jauh bet bet bet

Scatter Plot

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

Corrplot

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

Uji Korelasi antar Peubahs

cor.test(data$IPM, data$X1)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X1
## t = -2.613, df = 54, p-value = 0.0116
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.54952839 -0.07909634
## sample estimates:
##        cor 
## -0.3350301

cor.test(data$IPM, data$X2)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X2
## t = 17.48, df = 54, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8697444 0.9536370
## sample estimates:
##       cor 
## 0.9218567

cor.test(data$IPM, data$X3)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X3
## t = 12.516, df = 54, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7751966 0.9172900
## sample estimates:
##       cor 
## 0.8623505

cor.test(data$IPM, data$X4)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X4
## t = 5.8668, df = 54, p-value = 2.786e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4318579 0.7618552
## sample estimates:
##       cor 
## 0.6239211

cor.test(data$IPM, data$X5)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X5
## t = 11.436, df = 54, p-value = 4.795e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7426327 0.9041982
## sample estimates:
##       cor 
## 0.8412825

cor.test(data$IPM, data$X6)

## 
##  Pearson's product-moment correlation
## 
## data:  data$IPM and data$X6
## t = 4.7212, df = 54, p-value = 1.706e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3236139 0.7034623
## sample estimates:
##       cor 
## 0.5405273

Aman semua

Sebaran Spasial data

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

# Pastikan mapauas adalah objek sf
mapauas <- st_as_sf(mapauas)

peta <- st_make_valid(mapauas)

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

# Plot menggunakan ggplot2
ggplot(mapauas) +
  geom_sf(aes(fill = IPM)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Peta Sebaran Spasial Peubah Respon Y") +
  theme_minimal()

IPM yang tinggi cuman berpusat di Ibu kota- Ibu kotanya doang

Pembentukan Matrix Pembobot Spasial

Jadi ini keknya kan ada beberapa daerah yang terpisah pulau jadi ga cocok si Queen, tapi di coba dulu aje

Queen

# Memvalidasi geometri dalam objek sf
peta <- st_make_valid(mapauas)
# Konversi objek sf ke objek sp
sp.peta <- as(peta, "Spatial")

# Mengakses slot polygons dari objek sp
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 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

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

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.queen  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 2.5124, p-value = 0.005995
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.21143585       -0.01851852        0.00837704

Rook

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 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

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

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.rook  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 2.5124, p-value = 0.005995
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.21143585       -0.01851852        0.00837704

KNN

KNN 5

centroids <- st_centroid(mapauas)

## Warning: st_centroid assumes attributes are constant over geometries

longlat <- st_coordinates(centroids)



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

MI.knn5 <- moran(data$IPM,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mty = moran.test(data$IPM,W.knn.s,randomisation = T,alternative = "two.sided")
mty

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 3.4085, p-value = 0.0006533
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.238567510      -0.018181818       0.005674147

KNN 3

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

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

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 2.6237, p-value = 0.00435
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.239521273      -0.018181818       0.009647704

RDW

# Set zero.policy, false

W.rdw <-dnearneigh(longlat,0,85,longlat=TRUE)

## Warning in dnearneigh(longlat, 0, 85, longlat = TRUE): neighbour object has 19
## sub-graphs

W.rdw.s <- nb2listw(W.rdw,style='W',zero.policy=TRUE)

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

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.rdw.s  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.4139, p-value = 0.07869
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.17296140       -0.02500000        0.01960238

IDW

# Tambahkan koordinat long dan lat ke objek peta
peta$long <- longlat[, 1]  # Kolom X (long)
peta$lat  <- longlat[, 2]  # Kolom Y (lat

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


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$IPM,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(data$IPM,W.idw.s,randomisation = T,alternative = "greater")
mt3

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 3.5714, p-value = 0.0001775
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0668133643     -0.0181818182      0.0005663822

IDW Alpha 2

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$IPM,W.idw2.s,n=length(W.idw2.s$neighbours),S0=Szero(W.idw2.s))
mt4 = moran.test(data$IPM,W.idw2.s,randomisation = T,alternative = "greater")
mt4

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 2.8879, p-value = 0.001939
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.152435242      -0.018181818       0.003490476

EDW

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 #row-normalized
W.edw.s = mat2listw(W.edw.sd,style='W')

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

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 4.1778, p-value = 1.472e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.159123798      -0.018181818       0.001801142

EDW Alpha 2

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 #row-normalized
W.edw2.s = mat2listw(W.edw2.sd,style='W')
MI.edw2 <- moran(data$IPM,W.edw2.s,n=length(W.edw2.s$neighbours),S0=Szero(W.edw2.s))
mt6 = moran.test(data$IPM,W.edw2.s,randomisation = T,alternative = "greater")
mt6

## 
##  Moran I test under randomisation
## 
## data:  data$IPM  
## weights: W.edw2.s    
## 
## Moran I statistic standard deviate = 3.8071, p-value = 7.031e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.231809999      -0.018181818       0.004311909

Pemilihan

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.knn5$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.23856751 4.349522e-03
## 2      Radial Distance Weight (dmax=60)   0.23623996 7.869206e-02
## 3      Invers Distance Weight (alpha=1)   0.06681336 1.775333e-04
## 4      Invers Distance Weight (alpha=2)   0.15243524 1.939198e-03
## 5 Exponential Distance Weight (alpha=1)   0.15912380 1.471656e-05
## 6 Exponential Distance Weight (alpha=2)   0.23181000 7.031113e-05
## 7                       Rook Contiguity   0.21143585 5.994908e-03
## 8                      Queen Contiguity   0.21143585 5.994908e-03

KNN K= 5

Regresi OLS

ols<- lm(IPM~X1+X2+X3+X4+X5+X6,data=data)
summary(ols)

## 
## Call:
## lm(formula = IPM ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.27855 -0.24413  0.00432  0.35345  1.18372 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.519e+01  8.837e+00  -2.851 0.006361 ** 
## X1           1.582e-01  4.387e-02   3.606 0.000728 ***
## X2           1.562e+00  1.215e-01  12.859  < 2e-16 ***
## X3           7.366e-04  6.013e-05  12.249  < 2e-16 ***
## X4           8.384e-01  1.237e-01   6.776 1.47e-08 ***
## X5           9.541e-01  1.590e-01   6.002 2.32e-07 ***
## X6           6.770e-02  5.961e-02   1.136 0.261570    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5882 on 49 degrees of freedom
## Multiple R-squared:  0.9829, Adjusted R-squared:  0.9808 
## F-statistic: 469.9 on 6 and 49 DF,  p-value: < 2.2e-16

Peubah X1 tidak selaras dengan Corrplot dimana di Corrplot hubunganya negatif namun disini koefisiienya positif, sedangkan untuk peubah X6 tidak dimasukkan karena tidak signifikan. Akhirnya dipilihlah 3 peubah yaitu X2 X3 X4, karena X5 adalah Harapan Lama Sekolah yang dinilai kurang relevan dalam pembentukan IPM

Regresi

ols<- lm(IPM~X2+X3+X4,data=data)
summary(ols)

## 
## Call:
## lm(formula = IPM ~ X2 + X3 + X4, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4764 -0.5171 -0.1458  0.6656  2.5023 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.946e+01  1.153e+01  -2.554   0.0136 *  
## X2           2.014e+00  1.409e-01  14.297  < 2e-16 ***
## X3           7.140e-04  7.310e-05   9.768 2.33e-13 ***
## X4           1.031e+00  1.621e-01   6.360 5.18e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8376 on 52 degrees of freedom
## Multiple R-squared:  0.9632, Adjusted R-squared:  0.9611 
## F-statistic: 454.2 on 3 and 52 DF,  p-value: < 2.2e-16

Uji Asumsi OLS

Normalitas

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

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

hist(err.ols)

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

Homogenitas

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 0.19308, df = 3, p-value = 0.9787

Multikol

car::vif(ols)

##       X2       X3       X4 
## 2.181354 2.098171 1.321495

Autkorelasi Spasial

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 = IPM ~ X2 + X3 + X4, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 2.6344, p-value = 0.008428
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.206566086     -0.040179683      0.008772529

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 = 2.307, p-value = 0.02105
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.206566086      -0.018181818       0.009490393

Tolak H-0 Artinya ada autkorelasi spasial disini, Namun Normalitas masih terlanggar, maka dari itu dilakukan transformasi pada Peubah X3 ( Pengeluaran/ pendapatan ya gw lupa hehe) , agar skalanya sama

Transformasi Data

# Trnasformasi Logaritma X3
data$X3 <- log(data$X3)

ols<- lm(IPM~X2+X3+X4,data=data)
summary(ols)

## 
## Call:
## lm(formula = IPM ~ X2 + X3 + X4, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6605 -0.5562 -0.1007  0.6075  2.5424 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -100.4600    12.8177  -7.838 2.30e-10 ***
## X2             2.0747     0.1397  14.848  < 2e-16 ***
## X3             7.8412     0.8233   9.524 5.45e-13 ***
## X4             1.1056     0.1636   6.757 1.21e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8512 on 52 degrees of freedom
## Multiple R-squared:  0.962,  Adjusted R-squared:  0.9598 
## F-statistic: 439.2 on 3 and 52 DF,  p-value: < 2.2e-16

vif(ols)

##       X2       X3       X4 
## 2.077206 1.922150 1.304330

#normalitas dengan anderson darling
err.ols <- residuals(ols)
ad.test(err.ols)

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

#uji durbin watson
lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 0.30928, df = 3, p-value = 0.9583

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 = IPM ~ X2 + X3 + X4, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 3.1773, p-value = 0.001487
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.257209867     -0.040542663      0.008782145

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 = 2.838, p-value = 0.00454
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.257209867      -0.018181818       0.009416413

Aman ye jadinye, tinggal nyelasain autkorelasi spasialnya aja

Pemodelan Regspasi

ww= W.knn.s
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 = IPM ~ X2 + X3 + X4, data = data)
## test weights: listw
##  
##          statistic parameter  p.value   
## RSerr    6.3510639         1 0.011731 * 
## adjRSerr 6.7725479         1 0.009257 **
## RSlag    0.0068174         1 0.934196   
## adjRSlag 0.4283014         1 0.512824   
## SARMA    6.7793653         2 0.033719 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Signifikannya di SEM sama SARMA jadi dimodelin perbandingannya pake SEM Ddan sARMA

SEM

sem <- errorsarlm(IPM~X2+X3+X4, data=mapauas, listw=ww, zero.policy = TRUE)
summary(sem, Nagelkerke=T)

## 
## Call:errorsarlm(formula = IPM ~ X2 + X3 + X4, data = mapauas, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.3359582 -0.4693798  0.0062517  0.3741338  2.2634004 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -1.4125e+01  1.1766e+01 -1.2005    0.2299
## X2           2.0752e+00  1.4526e-01 14.2867 < 2.2e-16
## X3           7.5436e-04  7.4893e-05 10.0725 < 2.2e-16
## X4           8.0911e-01  1.6635e-01  4.8638 1.151e-06
## 
## Lambda: 0.39674, LR test value: 5.2218, p-value: 0.022305
## Asymptotic standard error: 0.13781
##     z-value: 2.879, p-value: 0.0039899
## Wald statistic: 8.2884, p-value: 0.0039899
## 
## Log likelihood: -64.84737 for error model
## ML residual variance (sigma squared): 0.56728, (sigma: 0.75318)
## Nagelkerke pseudo-R-squared: 0.96651 
## Number of observations: 56 
## Number of parameters estimated: 6 
## AIC: 141.69, (AIC for lm: 144.92)

lambdanya aman, intercept semua aman, rsquared juga bagus

ASUMSI

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

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

bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.24257, df = 3, p-value = 0.9704

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

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.095995, p-value = 0.9235
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.027518369      -0.018181818       0.009459592

aman tertatasi

SARMA

SARMA <- sacsarlm(IPM ~ X2+X3+X4,data=data,ww, zero.policy = TRUE)
summary(SARMA)

## 
## Call:sacsarlm(formula = IPM ~ X2 + X3 + X4, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.266660 -0.388155 -0.028901  0.443556  2.013831 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -78.22890   13.75785 -5.6861 1.300e-08
## X2            2.13640    0.13981 15.2803 < 2.2e-16
## X3            8.74563    0.82038 10.6604 < 2.2e-16
## X4            0.72902    0.16494  4.4199 9.873e-06
## 
## Rho: -0.047743
## Asymptotic standard error: 0.04473
##     z-value: -1.0674, p-value: 0.28581
## Lambda: 0.53874
## Asymptotic standard error: 0.11918
##     z-value: 4.5202, p-value: 6.1776e-06
## 
## LR test value: 10.382, p-value: 0.0055677
## 
## Log likelihood: -63.17523 for sac model
## ML residual variance (sigma squared): 0.5106, (sigma: 0.71456)
## Number of observations: 56 
## Number of parameters estimated: 7 
## AIC: 140.35, (AIC for lm: 146.73)

rhonya ga sinfinikan, lambdanya ga singikan, p-val semuanya aman sie

ASUMSI

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

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

bptest.Sarlm(SARMA)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.20385, df = 3, p-value = 0.977

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

## 
##  Moran I test under randomisation
## 
## data:  err.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.24428, p-value = 0.807
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.041936901      -0.018181818       0.009456997

LAH KALO GINI SARMA DONG

Pemilihan

df2 <- data.frame(
  "Model" = c("OLS (Regresi Klasik)", "SEM", "SARMA"),
  "AIC" = c(AIC(ols), AIC(sem), AIC(SARMA)),
  "Normalitas" = c(ad.test(err.ols)$p.value, ad.test(err.sem)$p.value, ad.test(err.SARMA)$p.value),
  "Heterogenitas" = c(lmtest::bptest(ols)$p.value, bptest.Sarlm(sem)$p.value, bptest.Sarlm(SARMA)$p.value),
  "Autokorelasi Spasial" = c(
    moran.test(err.ols, ww, alternative = "two.sided")$p.value,
    moran.test(err.sem, ww, alternative = "two.sided")$p.value,
    moran.test(err.SARMA, ww, alternative = "two.sided")$p.value
  )
)

df2

##                  Model      AIC Normalitas Heterogenitas Autokorelasi.Spasial
## 1 OLS (Regresi Klasik) 146.7320  0.1261584     0.9582746           0.00454012
## 2                  SEM 141.6947  0.5240016     0.9704417           0.92352425
## 3                SARMA 140.3505  0.3019873     0.9769658           0.80701745

SARMA

SARMA <- sacsarlm(IPM ~ X2+X3+X4,data=data,ww, zero.policy = TRUE)
summary(SARMA)

## 
## Call:sacsarlm(formula = IPM ~ X2 + X3 + X4, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.266660 -0.388155 -0.028901  0.443556  2.013831 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -78.22890   13.75785 -5.6861 1.300e-08
## X2            2.13640    0.13981 15.2803 < 2.2e-16
## X3            8.74563    0.82038 10.6604 < 2.2e-16
## X4            0.72902    0.16494  4.4199 9.873e-06
## 
## Rho: -0.047743
## Asymptotic standard error: 0.04473
##     z-value: -1.0674, p-value: 0.28581
## Lambda: 0.53874
## Asymptotic standard error: 0.11918
##     z-value: 4.5202, p-value: 6.1776e-06
## 
## LR test value: 10.382, p-value: 0.0055677
## 
## Log likelihood: -63.17523 for sac model
## ML residual variance (sigma squared): 0.5106, (sigma: 0.71456)
## Number of observations: 56 
## Number of parameters estimated: 7 
## AIC: 140.35, (AIC for lm: 146.73)

Cuman perlu dicatat RHonya ga singinikfna.

Spill over

Im = impacts(SARMA, listw = ww)
Im

## Impact measures (sac, exact):
##       Direct    Indirect     Total
## X2 2.1376033 -0.09855014 2.0390532
## X3 8.7505414 -0.40342706 8.3471143
## X4 0.7294326 -0.03362910 0.6958035

# Efek Umpan Balik
sum <- summary(SARMA)
sum

## 
## Call:sacsarlm(formula = IPM ~ X2 + X3 + X4, data = data, listw = ww, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.266660 -0.388155 -0.028901  0.443556  2.013831 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -78.22890   13.75785 -5.6861 1.300e-08
## X2            2.13640    0.13981 15.2803 < 2.2e-16
## X3            8.74563    0.82038 10.6604 < 2.2e-16
## X4            0.72902    0.16494  4.4199 9.873e-06
## 
## Rho: -0.047743
## Asymptotic standard error: 0.04473
##     z-value: -1.0674, p-value: 0.28581
## Lambda: 0.53874
## Asymptotic standard error: 0.11918
##     z-value: 4.5202, p-value: 6.1776e-06
## 
## LR test value: 10.382, p-value: 0.0055677
## 
## Log likelihood: -63.17523 for sac model
## ML residual variance (sigma squared): 0.5106, (sigma: 0.71456)
## Number of observations: 56 
## Number of parameters estimated: 7 
## AIC: 140.35, (AIC for lm: 146.73)

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

##    Koefisien EfekLangsung   UmpanBalik
## X2 2.1364033    2.1376033 0.0012000004
## X3 8.7456290    8.7505414 0.0049123487
## X4 0.7290231    0.7294326 0.0004094864

Peubah Penjelas X2

Pengaruh langsung dari peubah X1 adalah sebesar 2, artinya Jika rata-rata X1 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah tersebut akan meningkat sebesar 2 sekian , jika peubah lainnya tetap.
Pengaruh tidak langsung dari peubah X1 bernilai -0.0012, artinya jika Jika rata-rata X1 di wilayah-i meningkat satu satuan, maka rata-rata Y di wilayah-j (wilayah yang berbeda) akan berkurang sebesar -0.0012, jika peubah lainnya tetap.
Pengaruh total dari peubah X1 bernilai 1,98, 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 1.98
Notes disini Direct sama Indirectnya berbeda , perlu dilogikain dalam interpretasinya sih.

Kesimpulan

Santai dulu gak sih temen-temen hehe.