Packages

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.2
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.3
## Loading required package: carData
## 
## 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.2
## 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
library(lwgeom)       # Untuk memperbaiki geometri
## Warning: package 'lwgeom' was built under R version 4.3.3
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.11.2, PROJ 9.3.1
## 
## Attaching package: 'lwgeom'
## The following object is masked from 'package:sf':
## 
##     st_perimeter
library(stringr)
library(tmap)
## Warning: package 'tmap' was built under R version 4.3.3
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

Import Data

Data yang digunakan adalah data provinsi di Indonesia. Peubah - peubah yang digunakan adalah Y = Persentase penduduk miskin X1 = Persentase Rumah Tangga Sumber Penerangan Utama dari Listrik X2 = Persentase Rumah Tangga Akses Sumber Air Minum Layak X3 = Persentase Rumah Tangga Akses Sanitasi Layak X4 = Persentase Rumah Tangga Akses Hunian Layak X5 = Persentase Unmeet Pelayanan Kesehatan

Data

data <- read_excel("D:/Pekan Analisis Data/data.xlsx")
data$PROVINSI <- str_to_title(data$PROVINSI)
head(data)
## # A tibble: 6 × 7
##   PROVINSI           PPM  SPUL ASAML   ASL   AHL    PK
##   <chr>            <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aceh             14.2   99.9  90.1  81.1  67.8  5.24
## 2 Sumatera Utara    7.99  99.7  92.9  85.7  73.5  4.62
## 3 Sumatera Barat    5.97  99.4  86.8  72.8  62.3  4.8 
## 4 Riau              6.67  99.5  92.2  86.3  74.8  4.82
## 5 Jambi             7.1   99.6  82.2  84.0  66.2  4.84
## 6 Sumatera Selatan 11.0   99.6  87.2  82.4  63.2  4.92
data$PROVINSI <- str_replace(data$PROVINSI, "^Dki Jakarta$", "DKI Jakarta")

Dimensi Data

dim(data)
## [1] 38  7

Data terdiri dari 38 provinsi di Indonesia

Peta Kab/Kota

Berikut ini adalah peta Indonesia:

peta <- st_read("D:/Pekan Analisis Data/[LapakGIS.com] Batas Wilayah Provinsi 2024")
## Reading layer `LapakGIS_Batas_Provinsi_2024' from data source 
##   `D:\Pekan Analisis Data\[LapakGIS.com] Batas Wilayah Provinsi 2024' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY, XYZ
## Bounding box:  xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

Eksplorasi Data

Cek Missing Value

data[!complete.cases(data),]
## # A tibble: 0 × 7
## # ℹ 7 variables: PROVINSI <chr>, PPM <dbl>, SPUL <dbl>, ASAML <dbl>, ASL <dbl>,
## #   AHL <dbl>, PK <dbl>

Hasil pengecekan missing value menunjukkan bahwa tidak terdapat missing value pada data.

Ringkasan Staristik Deskriptif

summary(data)
##    PROVINSI              PPM             SPUL            ASAML      
##  Length:38          Min.   : 4.00   Min.   : 56.08   Min.   :30.64  
##  Class :character   1st Qu.: 6.32   1st Qu.: 99.12   1st Qu.:82.22  
##  Mode  :character   Median :10.13   Median : 99.72   Median :89.02  
##                     Mean   :11.15   Mean   : 97.01   Mean   :87.00  
##                     3rd Qu.:14.06   3rd Qu.: 99.92   3rd Qu.:94.65  
##                     Max.   :32.97   Max.   :100.00   Max.   :99.96  
##       ASL             AHL              PK       
##  Min.   :12.61   Min.   : 4.44   Min.   :1.390  
##  1st Qu.:79.38   1st Qu.:58.23   1st Qu.:4.112  
##  Median :83.36   Median :65.32   Median :4.830  
##  Mean   :81.14   Mean   :61.66   Mean   :4.870  
##  3rd Qu.:87.18   3rd Qu.:71.83   3rd Qu.:5.657  
##  Max.   :96.83   Max.   :86.68   Max.   :8.240
# Menunjukkan kab/kota dengan Y/PPM tertinggi
wilayah_max <- data %>%
  arrange(desc(data$PPM)) %>%
  slice(1)

print(paste("Provinsi dengan Y tertinggi adalah", wilayah_max$`PROVINSI`))
## [1] "Provinsi dengan Y tertinggi adalah Papua Pegunungan"
# Menunjukkan kab/kota dengan Y/PPM terendah
wilayah_min <- data %>%
  arrange(data$PPM) %>%
  slice(1)

print(paste("Provinsi dengan Y terendah adalah", wilayah_min$`PROVINSI`))
## [1] "Provinsi dengan Y terendah adalah Bali"

Berdasarkan output di atas diketahui bahwa Y/PPM tertinggi di Indonesia adalah Papua Pegunungan sebesar 32.97, sedangkan Y/PPM terendah adalah Bali sebesar 4.00. Rata-rata Y/PPM Provinsi di Indonesia adalah 11.15.

Boxplot

par(mfrow=c(1,6))
boxplot(data$PPM,main = "Sebaran Y/PPM", col="dodgerblue")
boxplot(data$SPUL,main = "Sebaran X1/SPUL", col="dodgerblue")
boxplot(data$ASAML,main = "Sebaran X2/ASAML", col="dodgerblue")
boxplot(data$ASL,main = "Sebaran X3/ASL", col="dodgerblue")
boxplot(data$AHL,main = "Sebaran X4/AHL", col="dodgerblue")
boxplot(data$AHL,main = "Sebaran X5/PK", col="dodgerblue")

Berdasarkan boxplot di atas dapat diketahui bahwa pada PPM, SPUL, ASAML, ASL, AHL, dan PK terdapat pencilan.

Hubungan antara Peubah Respon dan Peubah Penjelas

Scatter plot ini untuk melihat pola hubungan Y dengan masing-masing peubah X1, X2, X3, X4, X5 serta pola hubungan antar peubah X.

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

Berdasarkan scatter plot di atas, hubungan antara peubah independent dengan dependent serta antar peubah independent kurang terlihat jelas polanya.

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

Berdasarkan scatter plot ini, diperoleh nilai korelasi antara PPM atau peubah dependent dengan peubah independent yaitu ada hubungan negatif. Lalu, hubungan antar peubah independent bernilai positif.

cor(data[,2:7])
##              PPM       SPUL      ASAML        ASL        AHL         PK
## PPM    1.0000000 -0.7497862 -0.5719840 -0.7847651 -0.5738265 -0.2364989
## SPUL  -0.7497862  1.0000000  0.5413916  0.8292341  0.6840485  0.5069203
## ASAML -0.5719840  0.5413916  1.0000000  0.8044668  0.7083255  0.2305118
## ASL   -0.7847651  0.8292341  0.8044668  1.0000000  0.7513640  0.3306884
## AHL   -0.5738265  0.6840485  0.7083255  0.7513640  1.0000000  0.2988937
## PK    -0.2364989  0.5069203  0.2305118  0.3306884  0.2988937  1.0000000

Sebaran Spasial Data

# Gabungkan Data
combined_data <- peta %>%
  left_join(data, by = c("WADMPR" = "PROVINSI"))
combined_data
## Simple feature collection with 38 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY, XYZ
## Bounding box:  xmin: 94.97191 ymin: -11.00762 xmax: 141.02 ymax: 6.076832
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
## First 10 features:
##    KDPPUM                    WADMPR                                METADATA
## 1      11                      Aceh TASWIL1000020230928_DATA_BATAS_PROVINSI
## 2      12            Sumatera Utara TASWIL1000020230928_DATA_BATAS_PROVINSI
## 3      13            Sumatera Barat TASWIL1000020230928_DATA_BATAS_PROVINSI
## 4      14                      Riau TASWIL1000020230928_DATA_BATAS_PROVINSI
## 5      15                     Jambi TASWIL1000020230928_DATA_BATAS_PROVINSI
## 6      16          Sumatera Selatan TASWIL1000020230928_DATA_BATAS_PROVINSI
## 7      17                  Bengkulu TASWIL1000020230928_DATA_BATAS_PROVINSI
## 8      18                   Lampung TASWIL1000020230928_DATA_BATAS_PROVINSI
## 9      19 Kepulauan Bangka Belitung TASWIL1000020230928_DATA_BATAS_PROVINSI
## 10     21            Kepulauan Riau TASWIL1000020230928_DATA_BATAS_PROVINSI
##             UPDATED   PPM  SPUL ASAML   ASL   AHL   PK
## 1  Lapak GIS - 2024 14.23 99.86 90.08 81.10 67.76 5.24
## 2  Lapak GIS - 2024  7.99 99.73 92.94 85.73 73.47 4.62
## 3  Lapak GIS - 2024  5.97 99.39 86.85 72.82 62.29 4.80
## 4  Lapak GIS - 2024  6.67 99.52 92.24 86.32 74.80 4.82
## 5  Lapak GIS - 2024  7.10 99.58 82.16 83.97 66.21 4.84
## 6  Lapak GIS - 2024 10.97 99.61 87.23 82.36 63.21 4.92
## 7  Lapak GIS - 2024 13.56 99.88 72.10 83.01 56.52 6.35
## 8  Lapak GIS - 2024 10.69 99.94 84.51 85.44 65.96 5.19
## 9  Lapak GIS - 2024  4.55 99.87 82.97 94.16 30.72 4.78
## 10 Lapak GIS - 2024  5.37 99.91 94.71 91.23 57.19 4.61
##                          geometry
## 1  MULTIPOLYGON (((97.38631 1....
## 2  MULTIPOLYGON (((98.51629 -0...
## 3  MULTIPOLYGON (((100.7819 -3...
## 4  MULTIPOLYGON (((103.5003 -0...
## 5  MULTIPOLYGON (((104.4988 -1...
## 6  MULTIPOLYGON (((106.2155 -3...
## 7  MULTIPOLYGON (((102.385 -5....
## 8  MULTIPOLYGON (((105.4175 -6...
## 9  MULTIPOLYGON (((108.0467 -3...
## 10 MULTIPOLYGON (((105.2663 -1...

Sebaran Spasial Y/PPM

tmap_options(check.and.fix = TRUE)

# Mengatur mode peta
tmap_mode("plot")  # Ganti menjadi "view" jika Anda ingin peta interaktif
## tmap mode set to plotting
# Membuat peta dengan jumlah PPM dengan label nama provinsi
tm_shape(combined_data) + 
  tm_polygons("PPM", style = "quantile", n = 4, 
              title = "Jumlah PPM Indonesia") + 
  tm_text("WADMPR", size = 0.4, col = "black") +  # Tambahkan nama provinsi
  tm_layout(title = "Peta Jumlah PPM Indonesia", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol dengan wilayah tetangganya pada peubah PPM/Y di Indonesia. Hal ini tampak dari gradasi warna yang cenderung mengumpul. Artinya terdapat kemiripan antar kabupaten/kota terhadap sebaran PPM/Y sehingga terindikasi adanya autokorelasi spasial positif.

Matriks Pembobot Spasial

Matriks pembobot spasial dapat berdasarkan aspek ketetanggaan dan aspek jarak. Matriks ketetanggaan berdasarkan aspek jarak digunakan k-nearest neighbor (KNN), radial distance weight (RDW), inverse distance weight (IDW), dan exponential distance weight (IDW).

# Memvalidasi geometri dalam objek sf
peta <- st_make_valid(peta)

# Konversi objek sf ke objek sp
sp.peta <- as(peta, "Spatial")
sp.polygons <- SpatialPolygons(sp.peta@polygons)

Matriks Pembobot Berdasarkan Jarak

# Menghitung centroid dari setiap unit spasial
centroids <- st_centroid(peta)
## Warning: st_centroid assumes attributes are constant over geometries
# Ekstrak koordinat centroid
longlat <- st_coordinates(centroids)

# 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)

K-Nearest Neighbor

Matriks bobot k-nearest neighbor menentukan tetangga berdasarkan sejumlah tetangga terdekat (k) yang dipilih untuk setiap titik atau wilayah.

Karakteristik:

Untuk setiap titik/wilayah, tetangga yang dipilih adalah titik/wilayah yang paling dekat (berdasarkan jarak Euclidean atau jarak lainnya).

Nilai matriks bobot diberikan 1 jika suatu titik berada dalam k tetangga terdekat dari titik lain, dan 0 jika tidak.

Contoh: Misalkan k = 3, artinya setiap titik akan memiliki 3 tetangga terdekat. Matriks bobot akan memiliki 1 di elemen yang mewakili tetangga terdekat, dan 0 di elemen lainnya.

Digunakan 3 tetangga terdekat (k=3).

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

MI.knn <- moran(data$PPM,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(data$PPM,W.knn.s,randomisation = T,alternative = "greater")
mt1
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 6.1727, p-value = 3.356e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.67391466       -0.02702703        0.01289460

Inverse Distance Weight

Matriks pembobot Inverse Distance Weight (IDW) menentukan bobot berdasarkan kebalikan dari jarak antar titik atau wilayah. Semakin dekat dua titik, semakin besar bobot yang diberikan.

Ditentukan nilai alpha=1 dan alpha=2

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$PPM,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt2 = moran.test(data$PPM,W.idw.s,randomisation = T,alternative = "greater")
mt2
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 9.2903, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.299481359      -0.027027027       0.001235166
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$PPM,W.idw2.s,n=length(W.idw2.s$neighbours),S0=Szero(W.idw2.s))
mt3 = moran.test(data$PPM,W.idw2.s,randomisation = T,alternative = "greater")
mt3
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 7.1522, p-value = 4.269e-13
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.574974547      -0.027027027       0.007084562

Exponential Distance Weight

Matriks exponential distance weight memberikan bobot yang menurun secara eksponensial seiring dengan peningkatan jarak antar unit spasial. Ini mirip dengan inverse distance weight, tetapi penurunan bobot terhadap jarak lebih tajam.

Ditentukan alpha=1 dan alpha=2

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$PPM,W.edw.s,n=length(W.edw.s$neighbours),S0=Szero(W.edw.s))
mt4 = moran.test(data$PPM,W.edw.s,randomisation = T,alternative = "greater")
mt4
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 5.8323, p-value = 2.733e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.70840251       -0.02702703        0.01590003
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$PPM,W.edw2.s,n=length(W.edw2.s$neighbours),S0=Szero(W.edw2.s))
mt5 = moran.test(data$PPM,W.edw2.s,randomisation = T,alternative = "greater")
mt5
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.edw2.s    
## 
## Moran I statistic standard deviate = 4.8279, p-value = 6.898e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.73618288       -0.02702703        0.02499006

Pemilihan Matriks Pembobot Terbaik

Matriks pembobot spasial yang dipilih yaitu matriks yang menghasilkan nilai autokorelasi spasial yang tertinggi. Pemilihan matriks bobot juga dapat didasarkan pada matriks bobot dengan nilai indeks moran yang signifikan pada taraf 5%.

MatriksBobot <- c("K-Nearest Neighbor (k=3)",
                  "Invers Distance Weight (alpha=1)",
                  "Invers Distance Weight (alpha=2)",
                  "Exponential Distance Weight (alpha=1)",
                  "Exponential Distance Weight (alpha=2)")
IndeksMoran <- c(MI.knn$I,MI.idw$I,MI.idw2$I,MI.edw$I,MI.edw2$I)
pv = c(mt1$p.value,mt2$p.value,mt3$p.value,mt4$p.value,mt5$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=3)    0.6739147 3.355824e-10
## 2      Invers Distance Weight (alpha=1)    0.2994814 7.689766e-21
## 3      Invers Distance Weight (alpha=2)    0.5749745 4.269154e-13
## 4 Exponential Distance Weight (alpha=1)    0.7084025 2.732938e-09
## 5 Exponential Distance Weight (alpha=2)    0.7361829 6.898209e-07
Matriks[Matriks$`p-value`< 0.05,]
##                           Matriks Bobot Indeks Moran      p-value
## 1              K-Nearest Neighbor (k=3)    0.6739147 3.355824e-10
## 2      Invers Distance Weight (alpha=1)    0.2994814 7.689766e-21
## 3      Invers Distance Weight (alpha=2)    0.5749745 4.269154e-13
## 4 Exponential Distance Weight (alpha=1)    0.7084025 2.732938e-09
## 5 Exponential Distance Weight (alpha=2)    0.7361829 6.898209e-07

Berdasakan hasil di atas, Matriks bobot yang memiliki nilai Indeks Moran yang signifikan pada taraf 5% adalah k-nearest neighboor (k=3), invers distance weight (alpha = 1 dan alpha = 2), exponential distance weight (alpha = 1 dan alpha = 2).

Matriks pembobot spasial yang memiliki Indeks Moran tertinggi adalah EDW alpha = 2 sehingga matriks tersebut dapat dipilih menjadi matriks pembobot spasial terbaik diantara matriks pembobot spasial lainnya. Selanjutnya matriks pembobot yang digunakan dalam analisis adalah queen contiguity.

Model Regresi Klasik (OLS)

Pemodelan Ordinary Least Square dan uji asumsi multikolinearitas, autokorelasi spasial, kehomogenan ragam, kenormalan sisaan.

Matriks Bobot yang digunakan adalah EDW alpha = 2.

Metode OLS

ols <- lm(PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data)
summary(ols)
## 
## Call:
## lm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3436 -2.5597 -0.9662  1.8611 10.6258 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 58.242957  13.315875   4.374 0.000121 ***
## SPUL        -0.339116   0.182592  -1.857 0.072503 .  
## ASAML        0.002682   0.115109   0.023 0.981554    
## ASL         -0.239899   0.128871  -1.862 0.071874 .  
## AHL          0.033348   0.070465   0.473 0.639248    
## PK           0.611394   0.557528   1.097 0.280993    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.211 on 32 degrees of freedom
## Multiple R-squared:  0.6631, Adjusted R-squared:  0.6105 
## F-statistic:  12.6 on 5 and 32 DF,  p-value: 8.368e-07
print(paste("Nilai AIC OLS", AIC(ols)))
## [1] "Nilai AIC OLS 224.566279993226"

Model Diagnostics

Uji Normalitas Sisaan

H0 : sisaan model menyebar normal H1 : sisaan model tidak menyebar normal

Tolak H0 jika p-value < 0.05

err.ols <- residuals(ols)
ad.test(err.ols)
## 
##  Anderson-Darling normality test
## 
## data:  err.ols
## A = 0.46803, p-value = 0.2362

Karena p-value > 5%, maka tak tolak H0. Artinya sisaan menyebar normal.

hist(err.ols)

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

Berdasarkan plot histogram dan normal Q-Q plot diatas terlihat bahwa sisaan menyebar normal.

Uji Kehomogenan Ragam

H0 : Ragam Sisaan Homogen H1 : Ragam Sisaan Tidak Homogen

lmtest::bptest(ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 2.5996, df = 5, p-value = 0.7614

Karena p-value > 0.05, maka gagal tolak H0. Artinya ragam sisaan homogen

Uji Multikolinearitas

Pada regresi linear berganda perlu diperiksa persyaratan tidak adanya multikolinearitas antar peubah penjelas agar koefisien regresi yang diperoleh dapat dikatakan valid.

Pemeriksaan multikolinearitas dapat dilakukan dengan menggunakan nilai variance inflation factor (VIF). Multikolinearitas terjadi jika nilai VIF > 10.

car::vif(ols)
##     SPUL    ASAML      ASL      AHL       PK 
## 5.105869 3.863331 7.872506 2.714227 1.418314

Berdasarkan perhitungan, diperoleh nilai VIF < 10 yang menunjukkan tidak terjadi multikolinearitas antar peubah penjelas.

Uji Autokorelasi

H0 : Antargalat tidak memiliki autokorelasi H1 : Antargalat berkorelasi

dwtest(ols)
## 
##  Durbin-Watson test
## 
## data:  ols
## DW = 1.3879, p-value = 0.01231
## alternative hypothesis: true autocorrelation is greater than 0

Karena p-value < 0.05, maka tolak H0. Artinya, antargalat berkorelasi.

Pengujian Efek Spasial

Efek Dependensi Spasial

Untuk menentukan model dependensi spasial yang sesuai dengan data, terlebih dahulu dilakukan uji dependensi spasial. Secara umum ada dua macam uji dependensi spasia, yaitu Indeks Moran dan Uji Pengganda Lagrange (Lagrange Multiplier Test).

Uji Indeks Moran

Pengujian autokorelasi spasial menggunakan Indeks Moran yang memerlukan matriks pembobot spasial. Dalam kasus ini, digunakan matriks bobot terbaik yang telah diperoleh sebelumnya, yaitu radial distance weight.

H0 : Tidak ada autokorelasi spasial

H1 : Ada autokorelasi spasial

Uji Indeks Moran pada Galat Model Regresi

ww = W.edw2.s
lm.morantest(ols, listw=ww, alternative="two.sided")
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 3.2605, p-value = 0.001112
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.46650254      -0.06083937       0.02615914
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 = 3.0237, p-value = 0.002497
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.46650254       -0.02702703        0.02664139

Uji Indeks Moran pada peubah Y/PPM

my <- moran.test(data$PPM, ww,randomisation=T, 
           alternative="two.sided")
my
## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: ww    
## 
## Moran I statistic standard deviate = 4.8279, p-value = 1.38e-06
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.73618288       -0.02702703        0.02499006
moran.plot(data$PPM, ww, labels=data$`PROVINSI`)

Uji Indeks Moran pada peubah X1/SPUL

mx1 <- moran.test(data$SPUL, ww,randomisation=T, 
           alternative="two.sided")
mx1
## 
##  Moran I test under randomisation
## 
## data:  data$SPUL  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.2416, p-value = 0.001188
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.38293830       -0.02702703        0.01599427

Uji Indeks Moran pada peubah X2/ASAML

mx2 <- moran.test(data$ASAML, ww,randomisation=T, 
           alternative="two.sided")
mx2
## 
##  Moran I test under randomisation
## 
## data:  data$ASAML  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.9812, p-value = 6.856e-05
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.49768545       -0.02702703        0.01737030

Uji Indeks Moran pada peubah X3/ASL

mx3 <- moran.test(data$ASL, ww,randomisation=T, 
           alternative="two.sided")
mx3
## 
##  Moran I test under randomisation
## 
## data:  data$ASL  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.2892, p-value = 0.001005
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.41871155       -0.02702703        0.01836497

Uji Indeks Moran pada peubah X4/AHL

mx4 <- moran.test(data$AHL, ww,randomisation=T, 
           alternative="two.sided")
mx4
## 
##  Moran I test under randomisation
## 
## data:  data$AHL  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.036, p-value = 0.002397
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.44697331       -0.02702703        0.02437491

Uji Indeks Moran pada peubah X4/AHL

mx5 <- moran.test(data$PK, ww,randomisation=T, 
           alternative="two.sided")
mx5
## 
##  Moran I test under randomisation
## 
## data:  data$PK  
## weights: ww    
## 
## Moran I statistic standard deviate = 1.7693, p-value = 0.07684
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.26136794       -0.02702703        0.02656852

Berikut ini adalah hasil uji autokorelasi pada peubah respon, peubah penjelas, dan sisaan regresi

Peubah = c("Y", "X1", "X2","X3","X4", "X5", "Sisaan")
Indeks_Moran = c(my$estimate[1], mx1$estimate[1], mx2$estimate[1],mx3$estimate[1],mx4$estimate[1],mx5$estimate[1], merr$estimate[1])
p_value = c(my$p.value, mx1$p.value, mx2$p.value,mx3$p.value,mx4$p.value,mx5$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.7361829 1.379642e-06
## 2     X1    0.3829383 1.188435e-03
## 3     X2    0.4976855 6.855831e-05
## 4     X3    0.4187116 1.004866e-03
## 5     X4    0.4469733 2.397093e-03
## 6     X5    0.2613679 7.684213e-02
## 7 Sisaan    0.4665025 2.497267e-03

Berdasarkan hasil pengujian Indeks Moran di atas menunjukkan bahwa nilai indeks Moran pada peubah Y/PPM, X1/SPUL, X2/ASAML, X3/ASL, X4/AHL, dan Sisaan model regresi memiliki nilai p-value < 0.05 maka Tolak H0 sehingga autokorelasi spasial pada peubah tersebut berpengaruh nyata. Hal ini mengindikasikan bahwa terdapat ketergantungan spasial pada peubah Y/PPM, X1/SPUL, X2/ASAML, X3/ASL, X4/AHL, dan Sisaan regresi klasik. Terlihat pula nilai indeks moran positif pada peubah Y/PPM, X1/SPUL, X2/ASAML, X3/ASL, X4/AHL, dan Sisaan regresi klasik yang menunjukkan adanya autokorelasi positif.

Sedangkan, nilai indeks moran pada peubah X5/PK model memiliki p-value > 0.05 maka gagal tolak H0 artinya tidak terdapat autokorelasi spasial peubah X5/PK model taraf nyata 5%.

Oleh karenanya Spatial Durbin Watson dapat digunakan, untuk mencari model yang lebih baik, kita dapat melakukan uji LM (lagrange multiplier) untuk mengidentifikasi model dependensi spasial yang dapat digunakan pada kasus ini.

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 = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data)
## test weights: listw
##  
##          statistic parameter  p.value   
## RSerr      7.64295         1 0.005699 **
## RSlag     10.56213         1 0.001154 **
## adjRSerr   0.62579         1 0.428903   
## adjRSlag   3.54497         1 0.059726 . 
## SARMA     11.18792         2 0.003720 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Output menunjukkan bahwa hasil uji LM SAR, LM SEM, dan GSM signifikan pada taraf 5%. Berdasarkan skema tersebut, maka kita dapat mencoba kandidat model SAR, SEM, dan SARMA.

Efek Heterogenitas Spasial

H0 : Ragam Sisaan Homogen H1 : Ragam Sisaan Tidak Homogen

lmtest::bptest(ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 2.5996, df = 5, p-value = 0.7614

Karena p-value > 0.05, maka Gagal Tolak H0. Artinya ragam sisaan homogen atau tidak terdapat efek heterogenitas spasial.

Model Regresi Spasial

Model SAR

sar <- lagsarlm(PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, listw = ww, zero.policy=TRUE)
summary(sar, Nagelkerke = T)
## 
## Call:lagsarlm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, 
##     listw = ww, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8.19234 -2.24638  0.08671  2.49028  5.42334 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 33.665165  10.953708  3.0734 0.002116
## SPUL        -0.160679   0.137436 -1.1691 0.242356
## ASAML        0.050126   0.085482  0.5864 0.557612
## ASL         -0.262666   0.094135 -2.7903 0.005266
## AHL          0.024305   0.051436  0.4725 0.636547
## PK           0.863813   0.406385  2.1256 0.033536
## 
## Rho: 0.38691, LR test value: 14.58, p-value: 0.00013432
## Asymptotic standard error: 0.089672
##     z-value: 4.3148, p-value: 1.5975e-05
## Wald statistic: 18.617, p-value: 1.5975e-05
## 
## Log likelihood: -97.99306 for lag model
## ML residual variance (sigma squared): 9.4112, (sigma: 3.0678)
## Nagelkerke pseudo-R-squared: 0.77047 
## Number of observations: 38 
## Number of parameters estimated: 8 
## AIC: 211.99, (AIC for lm: 224.57)
## LM test for residual autocorrelation
## test value: 3.3386, p-value: 0.067672

Output di atas memperlihatkan bahwa koefisien Rho pada model SAR signifikan, dengan nilai AIC sebesar 211.99. Selain itu, terlihat pula hasil uji autokorelasi pada sisaan model memperlihatkan nilai p-value sebesar 0.067672, artinya tidak terdapat autokorelasi pada sisaan.

Uji Asumsi Model SAR

Asumsi Kenormalan Sisaan

H0 : galat model menyebar normal

H1 : galat model tidak menyebar normal

err.sar<-residuals(sar)
ad.test(err.sar)
## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.19144, p-value = 0.8907

Karena p-value > 5%, maka Gagal tolak H0. Artinya sisaan menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sar)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 4.0205, df = 5, p-value = 0.5465

Karena p-value > 5%, maka gagal tolak H0. Artinya ragam sisaan homogen.

Uji Kebebasan Sisaan

H0: Tidak ada autokorelasi spasial

H1: ada autokorelasi spasial

moran.test(err.sar, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: ww    
## 
## Moran I statistic standard deviate = 1.5418, p-value = 0.1231
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.22607997       -0.02702703        0.02694917

Karena p-value > 0.05, maka tolak H0 artinya tidak terdapat autokorelasi spasial pada sisaan model SAR.

Model SEM

sem <- errorsarlm(PPM ~ SPUL + ASAML + ASL + AHL + PK,data=data,listw=ww)
summary(sem)
## 
## Call:errorsarlm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, 
##     listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.58209 -1.95467 -0.31561  1.73582  6.35241 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 41.227373   9.685423  4.2566 2.075e-05
## SPUL        -0.192556   0.113687 -1.6937  0.090315
## ASAML       -0.051611   0.085123 -0.6063  0.544309
## ASL         -0.166330   0.079424 -2.0942  0.036241
## AHL          0.023492   0.046866  0.5013  0.616186
## PK           1.009599   0.343330  2.9406  0.003276
## 
## Lambda: 0.65674, LR test value: 18.328, p-value: 1.8593e-05
## Asymptotic standard error: 0.086396
##     z-value: 7.6015, p-value: 2.931e-14
## Wald statistic: 57.783, p-value: 2.931e-14
## 
## Log likelihood: -96.119 for error model
## ML residual variance (sigma squared): 7.0769, (sigma: 2.6603)
## Number of observations: 38 
## Number of parameters estimated: 8 
## AIC: 208.24, (AIC for lm: 224.57)

Output di atas menunjukkan bahwa koefisien Lambda tidak signifikan pada taraf nyata 5%. AIC model SEM adalah 208.24. Selanjutnya kita akan coba memeriksa sisaan model SEM ini.

Uji Asumsi Model SEM

Asumsi Kenormalan Sisaan

H0 : Sisaan model menyebar normal

H1 : Sisaaan model tidak menyebar normal

err.sem<-residuals(sem)
ad.test(err.sem)
## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.21306, p-value = 0.8418

Karena p-value > 0.05 maka Gagal tolak H0, artinya sisaan menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sem)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 4.1433, df = 5, p-value = 0.529

Karena p-value > 0.05 maka gagal tolak H0, artinya ragam sisaan homogen.

Asumsi Kebesaan Sisaan

H0 : Tidak Ada autokorelasi spasial

H1 : Ada autokorelasi spasial

moran.test(err.sem, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.51674, p-value = 0.6053
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.05794839       -0.02702703        0.02704200

Karena p-value > 0.05, maka gagal tolak H0 artinya tidak terdapat autokorelasi spasial pada sisaan model SEM.

Model SARMA/GSM

SARMA <- sacsarlm(PPM ~ SPUL + ASAML + ASL + AHL + PK,data=data,ww)
summary(SARMA)
## 
## Call:sacsarlm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, 
##     listw = ww)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -6.175643 -1.837702 -0.080191  1.824763  5.416249 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 40.417415  10.192280  3.9655 7.324e-05
## SPUL        -0.185379   0.121967 -1.5199   0.12853
## ASAML       -0.034577   0.087689 -0.3943   0.69335
## ASL         -0.211560   0.089535 -2.3629   0.01813
## AHL          0.025554   0.049278  0.5186   0.60406
## PK           1.040226   0.367441  2.8310   0.00464
## 
## Rho: 0.19019
## Asymptotic standard error: 0.13893
##     z-value: 1.3689, p-value: 0.17103
## Lambda: 0.50387
## Asymptotic standard error: 0.15341
##     z-value: 3.2845, p-value: 0.0010218
## 
## LR test value: 19.759, p-value: 5.1222e-05
## 
## Log likelihood: -95.4038 for sac model
## ML residual variance (sigma squared): 7.5847, (sigma: 2.754)
## Number of observations: 38 
## Number of parameters estimated: 9 
## AIC: 208.81, (AIC for lm: 224.57)

Output di atas memperlihatkan bahwa kedua koefisien dependensi spasial signifikan pada taraf nyata 5%, yaitu Rho dan Lambda. AIC model SARMA adalah sebesar 208.81.

Uji Asumsi Model SARMA

Asumsi Kenormalan Sisaan

H0 : Sisaan model menyebar normal

H1 : Sisaaan model tidak menyebar normal

err.SARMA<-residuals(SARMA)
ad.test(err.SARMA)
## 
##  Anderson-Darling normality test
## 
## data:  err.SARMA
## A = 0.12065, p-value = 0.9874

Karena p-value > 0.05 maka gagal tolak H0, artinya sisaan menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(SARMA)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 2.6771, df = 5, p-value = 0.7496

Karena p-value > 0.05 maka gagal tolak H0, artinya ragam sisaan homogen.

Asumsi Kebesaan Sisaan

H0 : Tidak Ada autokorelasi spasial

H1 : Ada autokorelasi spasial

moran.test(err.SARMA, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.42548, p-value = 0.6705
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.04318798       -0.02702703        0.02723369

Karena p-value > 0.05, maka gagal tolak H0 artinya tidak terdapat autokorelasi spasial pada sisaan model SARMA. Artinya model SARMA telah memenuhi asumsi kenormalan, kehomogenan ragam, dan kebebasan pada taraf nyata 5%.

Model SDM

sdm <- lagsarlm(PPM ~ SPUL + ASAML + ASL + AHL + PK, data=data, ww, zero.policy = TRUE, type="mixed");
summary(sdm)
## 
## Call:lagsarlm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, 
##     listw = ww, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.35527 -1.92268 -0.20279  1.43969  5.54898 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) 65.678609  21.289586  3.0850 0.002035
## SPUL        -0.344485   0.128628 -2.6782 0.007403
## ASAML       -0.078030   0.081975 -0.9519 0.341158
## ASL         -0.105710   0.087099 -1.2137 0.224871
## AHL          0.033839   0.044863  0.7543 0.450690
## PK           1.068596   0.353796  3.0204 0.002525
## lag.SPUL    -0.463348   0.253694 -1.8264 0.067789
## lag.ASAML    0.181652   0.118752  1.5297 0.126098
## lag.ASL      0.134338   0.146022  0.9200 0.357580
## lag.AHL      0.018551   0.059946  0.3095 0.756975
## lag.PK      -0.255646   0.450177 -0.5679 0.570117
## 
## Rho: 0.50014, LR test value: 11.484, p-value: 0.00070194
## Asymptotic standard error: 0.10714
##     z-value: 4.6683, p-value: 3.0374e-06
## Wald statistic: 21.793, p-value: 3.0374e-06
## 
## Log likelihood: -91.14651 for mixed model
## ML residual variance (sigma squared): 6.1848, (sigma: 2.4869)
## Number of observations: 38 
## Number of parameters estimated: 13 
## AIC: 208.29, (AIC for lm: 217.78)
## LM test for residual autocorrelation
## test value: 0.048771, p-value: 0.82522

Uji Asumsi Model SDM ### Asumsi Kenormalan Sisaan

H0 : Sisaan model menyebar normal

H1 : Sisaaan model tidak menyebar normal

err.sdm<-sdm$residuals 
ad.test(err.sdm) #anderson darling
## 
##  Anderson-Darling normality test
## 
## data:  err.sdm
## A = 0.25135, p-value = 0.7225

Karena p-value > 0.05 maka gagal tolak H0, artinya sisaan menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sdm)
## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 7.9998, df = 10, p-value = 0.6289

Karena p-value > 0.05 maka gagal tolak H0, artinya ragam sisaan homogen.

Asumsi Kebesaan Sisaan

H0 : Tidak Ada autokorelasi spasial

H1 : Ada autokorelasi spasial

moran.test(err.sdm, ww, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  err.sdm  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.25862, p-value = 0.7959
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.01568041       -0.02702703        0.02726915

Karena p-value > 0.05, maka gagal tolak H0 artinya tidak terdapat autokorelasi spasial pada sisaan model sdm Artinya model sdm telah memenuhi asumsi kenormalan, kehomogenan ragam, dan kebebasan pada taraf nyata 5%.

Goodness of Fits (Model Terbaik)

Model regresi terbaik ditentukan berdasarkan nilai AIC yang dihasilkan dari masing-masing model. Model terbaik akan menghasilkan nilai AIC yang kecil. Hasilnya adalah sebagai berikut:

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) 224.5663
## 2                  SAR 211.9861
## 3                  SEM 208.2380
## 4                SARMA 208.8076
## 5                  SDM 208.2930

Berdasarkan output diatas, model SEM adalah model yang terbaik berdasarkan nilai AIC-nya. Hal ini juga sejalan dengan Lambda yang signifikan (p-value 0.00). Kemudian, hasil uji asumsi sisaan juga menunjukkan bahwa model SEM telah memenuhi asumsi kenormalan, kehomogenan ragam, dan kebebasan.

Interpretasi dan Kesimpulan

sum <- summary(sem)
sum
## 
## Call:errorsarlm(formula = PPM ~ SPUL + ASAML + ASL + AHL + PK, data = data, 
##     listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.58209 -1.95467 -0.31561  1.73582  6.35241 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 41.227373   9.685423  4.2566 2.075e-05
## SPUL        -0.192556   0.113687 -1.6937  0.090315
## ASAML       -0.051611   0.085123 -0.6063  0.544309
## ASL         -0.166330   0.079424 -2.0942  0.036241
## AHL          0.023492   0.046866  0.5013  0.616186
## PK           1.009599   0.343330  2.9406  0.003276
## 
## Lambda: 0.65674, LR test value: 18.328, p-value: 1.8593e-05
## Asymptotic standard error: 0.086396
##     z-value: 7.6015, p-value: 2.931e-14
## Wald statistic: 57.783, p-value: 2.931e-14
## 
## Log likelihood: -96.119 for error model
## ML residual variance (sigma squared): 7.0769, (sigma: 2.6603)
## Number of observations: 38 
## Number of parameters estimated: 8 
## AIC: 208.24, (AIC for lm: 224.57)

Berikut ini adalah interpretasinya:

Peubah Penjelas X1/Sumber Penerangan Utama Listrik (SPUL)

Jika akses sumber penerangan utama listrik naik 1 persen maka persentase penduduk miskin turun 0.192556.

Peubah Penjelas X2/Akses Sumber Air Minum Layak

Jika akses sumber air minum layak naik 1 persen maka persentase penduduk miskin turun 0.192556.

Peubah Penjelas X3/Akses Sanitasi Layak (ASL)

Jika akses sanitasi layak naik 1 persen maka persentase penduduk miskin turun 0.166330.

Peubah Penjelas X4/Akses Hunian Layak

Jika akses hunian layak naik 1 persen maka persentase penduduk miskin naik 0.023492.

Peubah Penjelas X5/Pelayanan Kesehatan

Jika pelayanan kesehatan naik 1 persen maka persentase penduduk miskin naik 1.009599.

Lambda

Nilai lambda (0.65674) menunjukkan adanya autokorelasi spasial positif yang kuat dalam model, di mana pola residual di suatu wilayah memiliki hubungan yang erat dengan pola residual di wilayah sekitarnya. Dengan memasukkan lambda ke dalam model, SEM mampu menangkap pola ini, sehingga hasil estimasi menjadi lebih akurat dibandingkan model regresi konvensional.