Project Regresi Spasial

Packages

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

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

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

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.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.5.3

library(sf)           # Mengelola data spasial dengan standar Simple Features (sf)
library(readxl)       # Membaca file Excel (format .xls dan .xlsx)
library(openxlsx)     # Membuat, membaca, dan menulis file Excel (lebih cepat dari readxl)

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

library(corrplot)     # Membuat visualisasi korelasi dalam bentuk plot

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

## corrplot 0.95 loaded

library(DescTools)    # Kumpulan fungsi statistik deskriptif dan inferensial

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

## Loading required package: carData

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

## 
## 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.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)       # Fungsi untuk menguji model regresi (heteroskedastisitas, serial correlation)

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

## Loading required package: zoo

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

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

## 
## 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.5.3

Import Data

Data yang digunakan adalah data sosial ekonomi dari kab/kota di Provinsi Sumatera Utara. Peubah - peubah yang digunakan adalah Y = Persentase penduduk miskin X1 = Rata - rata lama sekolah X2 = Pengeluaran per kapita X3 = Umur harapan hidup X4 = Harapan lama sekolah

Data

data <- read_excel("C:\\Users\\ASUS\\Downloads\\Sumut.xlsx")
head(data)

## # A tibble: 6 × 8
##   `NAMA KABUPATEN`     PPM   RLS   PPK   UHH   HLS   lat  long
##   <chr>              <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Asahan              8.21  8.83 11795  73.4  12.6  2.82  99.6
## 2 Batu Bara          11.4   8.5  10933  72.6  13.1  3.17  99.5
## 3 Dairi               7.47  9.88 10969  74.1  13.3  2.87  98.3
## 4 Deli Serdang        3.44 10.3  12890  73.6  13.4  3.42  98.7
## 5 Humbang Hasundutan  8.69 10.0   8476  74.1  13.3  2.20  98.6
## 6 Karo                7.98 10.0  12779  74.2  13.2  3.11  98.3

Dimensi Data

dim(data)

## [1] 33  8

Data terdiri dari 33 kab/kota di Provinsi Sumatera Utara

Peta Kab/Kota

Berikut ini adalah peta Provinsi Sumatera Utara:

peta <- st_read("D:/Sem 5/RegSpas/Tugas/UAS/RBI_50K_2023_Sumatera Utara.x26272")

## Reading layer `RBI_50K_2023_Sumatera Utara' from data source 
##   `D:\Sem 5\RegSpas\Tugas\UAS\RBI_50K_2023_Sumatera Utara.x26272' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 97.05747 ymin: -0.6387974 xmax: 100.4345 ymax: 4.302547
## Geodetic CRS:  WGS 84

# Filter data untuk menghapus entri dengan NAMOBJ "Sumatera Utara"
peta <- subset(peta, NAMOBJ != "Sumatera Utara")

peta

## Simple feature collection with 33 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 97.05747 ymin: -0.6387974 xmax: 100.4345 ymax: 4.302547
## Geodetic CRS:  WGS 84
## First 10 features:
##                   NAMOBJ      FCODE REMARK                   METADATA SRS_ID
## 1                 Asahan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 2              Batu Bara BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 3                  Dairi BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 4           Deli Serdang BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 5     Humbang Hasundutan BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 6                   Karo BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 7            Kota Binjai BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 8      Kota Gunungsitoli BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 9             Kota Medan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 10 Kota Padang Sidempuan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
##    KDBBPS KDCBPS KDCPUM KDEBPS KDEPUM KDPBPS KDPKAB KDPPUM     LUASWH TIPADM
## 1    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.09     12 3737.82989      4
## 2    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.19     12  888.14237      4
## 3    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.11     12 2083.60430      4
## 4    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.07     12 2581.23236      4
## 5    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.16     12 2351.51422      4
## 6    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.06     12 2206.87627      4
## 7    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.75     12   93.77007      5
## 8    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.78     12  208.68440      5
## 9    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.71     12  279.29043      5
## 10   <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.77     12  159.29837      5
##    WADMKC WADMKD                WADMKK         WADMPR WIADKC           WIADKK
## 1    <NA>   <NA>                Asahan Sumatera Utara   <NA>             <NA>
## 2    <NA>   <NA>             Batu Bara Sumatera Utara   <NA>           Asahan
## 3    <NA>   <NA>                 Dairi Sumatera Utara   <NA>             <NA>
## 4    <NA>   <NA>          Deli Serdang Sumatera Utara   <NA>             <NA>
## 5    <NA>   <NA>    Humbang Hasundutan Sumatera Utara   <NA>   Tapanuli Utara
## 6    <NA>   <NA>                  Karo Sumatera Utara   <NA>             <NA>
## 7    <NA>   <NA>           Kota Binjai Sumatera Utara   <NA>             <NA>
## 8    <NA>   <NA>     Kota Gunungsitoli Sumatera Utara   <NA>             Nias
## 9    <NA>   <NA>            Kota Medan Sumatera Utara   <NA>             <NA>
## 10   <NA>   <NA> Kota Padang Sidempuan Sumatera Utara   <NA> Tapanuli Selatan
##    WIADPR WIADKD SHAPE_Leng  SHAPE_Area                       geometry
## 1    <NA>      0  3.7686895 0.304020059 MULTIPOLYGON (((99.73741 3....
## 2    <NA>      0  2.4191333 0.072265659 MULTIPOLYGON (((99.72279 3....
## 3    <NA>      0  2.3768296 0.169476903 MULTIPOLYGON (((98.57832 2....
## 4    <NA>      0  4.5635597 0.210079808 MULTIPOLYGON (((98.74528 3....
## 5    <NA>      0  2.3376370 0.191184267 MULTIPOLYGON (((98.84978 2....
## 6    <NA>      0  2.7342061 0.179547636 MULTIPOLYGON (((98.47139 3....
## 7    <NA>      0  0.6512809 0.007632627 MULTIPOLYGON (((98.49823 3....
## 8    <NA>      0  0.9416521 0.016957951 MULTIPOLYGON (((97.65328 1....
## 9    <NA>      0  1.7186495 0.022734196 MULTIPOLYGON (((98.74528 3....
## 10   <NA>      0  0.4740388 0.012945257 MULTIPOLYGON (((99.31969 1....

plot(peta, main = "Peta Provinsi Sumatera Utara")

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

plot(peta[1], main = "Peta Provinsi Sumatera Utara (Kab/Kota)")

Eksplorasi Data

Cek Missing Value

data[!complete.cases(data),]

## # A tibble: 0 × 8
## # ℹ 8 variables: NAMA KABUPATEN <chr>, PPM <dbl>, RLS <dbl>, PPK <dbl>,
## #   UHH <dbl>, HLS <dbl>, lat <dbl>, long <dbl>

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

Ringkasan Staristik Deskriptif

summary(data)

##  NAMA KABUPATEN          PPM              RLS              PPK       
##  Length:33          Min.   : 3.440   Min.   : 6.140   Min.   : 6382  
##  Class :character   1st Qu.: 7.870   1st Qu.: 8.840   1st Qu.: 9395  
##  Mode  :character   Median : 8.540   Median : 9.510   Median :11632  
##                     Mean   : 9.925   Mean   : 9.387   Mean   :10911  
##                     3rd Qu.:11.420   3rd Qu.:10.090   3rd Qu.:12115  
##                     Max.   :22.810   Max.   :11.620   Max.   :15674  
##       UHH             HLS             lat              long       
##  Min.   :71.52   Min.   :12.64   Min.   :0.7086   Min.   : 97.39  
##  1st Qu.:72.30   1st Qu.:13.25   1st Qu.:1.5759   1st Qu.: 98.31  
##  Median :73.65   Median :13.42   Median :2.3466   Median : 99.06  
##  Mean   :73.25   Mean   :13.48   Mean   :2.2944   Mean   : 98.89  
##  3rd Qu.:74.10   3rd Qu.:13.70   3rd Qu.:2.9782   3rd Qu.: 99.37  
##  Max.   :74.76   Max.   :14.78   Max.   :3.8654   Max.   :100.17

# Menunjukkan kab/kota dengan Y/PPM tertinggi
wilayah_max <- data %>%
  arrange(desc(data$PPM)) %>%
  slice(1)

print(paste("Kabupaten/kota dengan Y tertinggi adalah", wilayah_max$`NAMA KABUPATEN`))

## [1] "Kabupaten/kota dengan Y tertinggi adalah Nias Barat"

# Menunjukkan kab/kota dengan Y/PPM terendah
wilayah_min <- data %>%
  arrange(data$PPM) %>%
  slice(1)

print(paste("Kabupaten/kota dengan Y terendah adalah", wilayah_min$`NAMA KABUPATEN`))

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

Berdasarkan output di atas diketahui bahwa Y/PPM tertinggi di Provinsi Sumatera Utara adalah Nias Barat sebesar 22.81, sedangkan Y/PPM terendah adalah Deli Serdang sebesar 3.44. Rata-rata Y/PPM kab/kota di Provinsi Sumatera Utara adalah 9.925.

Boxplot

par(mfrow=c(1,5))
boxplot(data$PPM,main = "Sebaran Y/PPM", col="dodgerblue")
boxplot(data$RLS,main = "Sebaran X1/RLS", col="dodgerblue")
boxplot(data$PPK,main = "Sebaran X2/PPK", col="dodgerblue")
boxplot(data$UHH,main = "Sebaran X3/UHH", col="dodgerblue")
boxplot(data$HLS,main = "Sebaran X4/HLS", col="dodgerblue")

Berdasarkan boxplot di atas dapat diketahui bahwa pada PPM, RLS, dan HLS terdapat pencilan.

Hubungan antara Peubah Respon dan Peubah Penjelas

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

pairs(data[,2:6], 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:6]), 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:6])

##            PPM        RLS        PPK        UHH        HLS
## PPM  1.0000000 -0.7441833 -0.7214875 -0.3481624 -0.3068192
## RLS -0.7441833  1.0000000  0.7722856  0.6529525  0.5728017
## PPK -0.7214875  0.7722856  1.0000000  0.5691507  0.2973662
## UHH -0.3481624  0.6529525  0.5691507  1.0000000  0.2446364
## HLS -0.3068192  0.5728017  0.2973662  0.2446364  1.0000000

Sebaran Spasial Data

# Gabungkan Data
combined_data <- peta %>%
  left_join(data, by = c("NAMOBJ" = "NAMA KABUPATEN"))
combined_data

## Simple feature collection with 33 features and 32 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 97.05747 ymin: -0.6387974 xmax: 100.4345 ymax: 4.302547
## Geodetic CRS:  WGS 84
## First 10 features:
##                   NAMOBJ      FCODE REMARK                   METADATA SRS_ID
## 1                 Asahan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 2              Batu Bara BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 3                  Dairi BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 4           Deli Serdang BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 5     Humbang Hasundutan BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 6                   Karo BA03050040 <Null> TASWIL5000020230907KABKOTA   4326
## 7            Kota Binjai BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 8      Kota Gunungsitoli BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 9             Kota Medan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
## 10 Kota Padang Sidempuan BA03050040   <NA> TASWIL5000020230907KABKOTA   4326
##    KDBBPS KDCBPS KDCPUM KDEBPS KDEPUM KDPBPS KDPKAB KDPPUM     LUASWH TIPADM
## 1    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.09     12 3737.82989      4
## 2    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.19     12  888.14237      4
## 3    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.11     12 2083.60430      4
## 4    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.07     12 2581.23236      4
## 5    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.16     12 2351.51422      4
## 6    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.06     12 2206.87627      4
## 7    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.75     12   93.77007      5
## 8    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.78     12  208.68440      5
## 9    <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.71     12  279.29043      5
## 10   <NA>   <NA>   <NA>   <NA>   <NA>   <NA>  12.77     12  159.29837      5
##    WADMKC WADMKD                WADMKK         WADMPR WIADKC           WIADKK
## 1    <NA>   <NA>                Asahan Sumatera Utara   <NA>             <NA>
## 2    <NA>   <NA>             Batu Bara Sumatera Utara   <NA>           Asahan
## 3    <NA>   <NA>                 Dairi Sumatera Utara   <NA>             <NA>
## 4    <NA>   <NA>          Deli Serdang Sumatera Utara   <NA>             <NA>
## 5    <NA>   <NA>    Humbang Hasundutan Sumatera Utara   <NA>   Tapanuli Utara
## 6    <NA>   <NA>                  Karo Sumatera Utara   <NA>             <NA>
## 7    <NA>   <NA>           Kota Binjai Sumatera Utara   <NA>             <NA>
## 8    <NA>   <NA>     Kota Gunungsitoli Sumatera Utara   <NA>             Nias
## 9    <NA>   <NA>            Kota Medan Sumatera Utara   <NA>             <NA>
## 10   <NA>   <NA> Kota Padang Sidempuan Sumatera Utara   <NA> Tapanuli Selatan
##    WIADPR WIADKD SHAPE_Leng  SHAPE_Area   PPM   RLS   PPK   UHH   HLS    lat
## 1    <NA>      0  3.7686895 0.304020059  8.21  8.83 11795 73.39 12.64 2.8175
## 2    <NA>      0  2.4191333 0.072265659 11.38  8.50 10933 72.63 13.11 3.1741
## 3    <NA>      0  2.3768296 0.169476903  7.47  9.88 10969 74.13 13.32 2.8676
## 4    <NA>      0  4.5635597 0.210079808  3.44 10.28 12890 73.65 13.39 3.4202
## 5    <NA>      0  2.3376370 0.191184267  8.69 10.01  8476 74.07 13.32 2.1989
## 6    <NA>      0  2.7342061 0.179547636  7.98 10.03 12779 74.16 13.25 3.1053
## 7    <NA>      0  0.6512809 0.007632627  4.79 11.19 11567 74.18 14.17 3.6135
## 8    <NA>      0  0.9416521 0.016957951 14.78  8.65  8635 74.03 13.78 1.2805
## 9    <NA>      0  1.7186495 0.022734196  8.00 11.62 15674 74.76 14.78 3.5952
## 10   <NA>      0  0.4740388 0.012945257  6.85 11.12 11552 73.54 14.59 1.3722
##       long                       geometry
## 1  99.6341 MULTIPOLYGON (((99.73741 3....
## 2  99.5006 MULTIPOLYGON (((99.72279 3....
## 3  98.2651 MULTIPOLYGON (((98.57832 2....
## 4  98.7041 MULTIPOLYGON (((98.74528 3....
## 5  98.5721 MULTIPOLYGON (((98.84978 2....
## 6  98.2651 MULTIPOLYGON (((98.47139 3....
## 7  98.5025 MULTIPOLYGON (((98.49823 3....
## 8  97.6147 MULTIPOLYGON (((97.65328 1....
## 9  98.6722 MULTIPOLYGON (((98.74528 3....
## 10 99.2730 MULTIPOLYGON (((99.31969 1....

Sebaran Spasial Y/PPM

library(tmap)

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

## Registered S3 method overwritten by 'stars':
##   method                  from
##   st_interpolate_aw.stars sf

combined_data <- st_make_valid(combined_data)
tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.

# Membuat peta dengan jumlah PPM dengan label nama kabupaten/kota
tm_shape(combined_data) + 
  tm_polygons("PPM", style = "quantile", n = 4, 
              title = "Jumlah PPM Sumatera Utara") + 
  tm_text("NAMOBJ", size = 0.4, col = "black") +  # Tambahkan nama kabupaten/kota
  tm_layout(title = "Peta Jumlah PPM Sumatera Utara", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n' to 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol dengan wilayah tetangganya pada peubah PPM/Y di Provinsi Sumatera Utara. 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.

Sebaran Spasial X1/RLS

# Membuat peta dengan jumlah RLS dengan label nama kabupaten/kota
tm_shape(combined_data) + 
  tm_polygons("RLS", style = "quantile", n = 4, 
              title = "Jumlah RLS Sumatera Utara") + 
  tm_text("NAMOBJ", size = 0.4, col = "black") +  # Tambahkan nama kabupaten/kota
  tm_layout(title = "Peta Jumlah RLS Sumatera Utara", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Sebaran Spasial X2/PPK

# Membuat peta dengan jumlah PPK dengan label nama kabupaten/kota
tm_shape(combined_data) + 
  tm_polygons("PPK", style = "quantile", n = 4, 
              title = "Jumlah PPK Sumatera Utara") + 
  tm_text("NAMOBJ", size = 0.4, col = "black") +  # Tambahkan nama kabupaten/kota
  tm_layout(title = "Peta Jumlah PPK Sumatera Utara", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Sebaran Spasial X3/UHH

# Membuat peta dengan jumlah UHH dengan label nama kabupaten/kota
tm_shape(combined_data) + 
  tm_polygons("UHH", style = "quantile", n = 4, 
              title = "Jumlah UHH Sumatera Utara") + 
  tm_text("NAMOBJ", size = 0.4, col = "black") +  # Tambahkan nama kabupaten/kota
  tm_layout(title = "Peta Jumlah UHH Sumatera Utara", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Sebaran Spasial X4/HLS

# Membuat peta dengan jumlah HLS dengan label nama kabupaten/kota
tm_shape(combined_data) + 
  tm_polygons("HLS", style = "quantile", n = 4, 
              title = "Jumlah HLS Sumatera Utara") + 
  tm_text("NAMOBJ", size = 0.4, col = "black") +  # Tambahkan nama kabupaten/kota
  tm_layout(title = "Peta Jumlah HLS Sumatera Utara", 
            legend.outside = TRUE, 
            legend.title.size = 1, 
            legend.text.size = 0.8)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## map variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Matriks Pembobot Spasial

Matriks pembobot spasial dapat berdasarkan aspek ketetanggaan dan aspek jarak. Matriks ketetanggaan yang dicobakan dalam analisis ini yaitu rook contiguity dan queen contiguity sedangkan berdasarkan aspek jarak digunakan k-nearest neighbor (KNN), radial distance weight (RDW), inverse distance weight (IDW), dan exponential distance weight (IDW).

Matriks Pembobot Berdasarkan Ketetanggaan

Queen Contiguity

Matriks pembobot spasial queen contiguity lebih luas daripada rook contiguity karena wilayah-wilayah dianggap bertetangga jika mereka berbagi sisi atau titik sudut.

Karakteristik:

Dua poligon dianggap bertetangga jika mereka berbagi batas (baik sisi maupun titik sudut).

Jika dua wilayah memiliki titik atau sisi yang bersinggungan, maka elemen matriks pembobot diberi nilai 1.

# Memvalidasi geometri dalam objek sf
peta <- st_make_valid(peta)
# 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): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

queen

## Neighbour list object:
## Number of regions: 33 
## Number of nonzero links: 118 
## Percentage nonzero weights: 10.83563 
## Average number of links: 3.575758 
## 2 disjoint connected subgraphs

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

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.queen    
## 
## Moran I statistic standard deviate = 6.2127, p-value = 2.604e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.76773982       -0.03125000        0.01653951

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

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

Rook Contiguity

Matriks pembobot spasial rook contiguity didasarkan pada konsep bahwa dua wilayah dianggap bertetangga jika mereka berbagi sisi batas.

Karakteristik:

Dua poligon dianggap bertetangga jika mereka memiliki batas yang bersinggungan di sepanjang sisi yang sama (tidak hanya di titik sudut).

Jika dua wilayah bertetangga, nilai elemen matriks pembobot adalah 1; jika tidak, maka nilai elemen tersebut adalah 0.

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

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

rook

## Neighbour list object:
## Number of regions: 33 
## Number of nonzero links: 118 
## Percentage nonzero weights: 10.83563 
## Average number of links: 3.575758 
## 2 disjoint connected subgraphs

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

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.rook    
## 
## Moran I statistic standard deviate = 6.2127, p-value = 2.604e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.76773982       -0.03125000        0.01653951

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

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

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 5 tetangga terdekat (k=5).

# k = 5
W.knn<-knn2nb(knearneigh(longlat,k=5,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.9406, p-value = 1.952e-12
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.603119698      -0.031250000       0.008353978

Digunakan 3 tetangga terdekat (k=3).

# k = 3
W.knn.3 <-knn2nb(knearneigh(longlat,k=3,longlat=TRUE))

## Warning in knn2nb(knearneigh(longlat, k = 3, longlat = TRUE)): neighbour object
## has 2 sub-graphs

W.knn.3.s <- nb2listw(W.knn,style='W')

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

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 6.9406, p-value = 1.952e-12
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.603119698      -0.031250000       0.008353978

Radial Distance Weight

Matriks bobot radial distance mendefinisikan tetangga berdasarkan jarak geografis. Setiap titik atau wilayah dianggap sebagai tetangga jika jaraknya dari titik/wilayah pusat berada dalam jarak radius tertentu.

Karakteristik:

Tetangga dipilih berdasarkan jarak dari titik pusat hingga ke titik lain dalam radius tertentu (r). Jika titik lain berada dalam radius tersebut, maka nilai dalam matriks bobot adalah 1; jika tidak, maka 0.

Jarak dihitung berdasarkan metrik spasial seperti jarak Euclidean.

Contoh: Jika radius yang ditentukan adalah 10 km, maka setiap titik atau wilayah yang berada dalam jarak 10 km dari titik pusat akan dianggap sebagai tetangga, dan bobot matriks akan diberikan 1.

Ditentukan nilai ambang dmax adalah 105 km (d=105). Nilai ini menggambarkan nilai jarak maksimum untuk menentukan dependensi spasial anatar lokasi-i terhadap lokasi-j.

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

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

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

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

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.rdw.s    
## 
## Moran I statistic standard deviate = 9.5265, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.712965482      -0.031250000       0.006102799

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))
mt3 = moran.test(data$PPM,W.idw.s,randomisation = T,alternative = "greater")
mt3

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 7.2421, p-value = 2.209e-13
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.276568851      -0.031250000       0.001806597

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))
mt4 = moran.test(data$PPM,W.idw2.s,randomisation = T,alternative = "greater")
mt4

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 5.8423, p-value = 2.574e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.58469819       -0.03125000        0.01111513

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))
mt5 = moran.test(data$PPM,W.edw.s,randomisation = T,alternative = "greater")
mt5

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 10.703, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2656130663     -0.0312500000      0.0007692685

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))
mt6 = moran.test(data$PPM,W.edw2.s,randomisation = T,alternative = "greater")
mt6

## 
##  Moran I test under randomisation
## 
## data:  data$PPM  
## weights: W.edw2.s    
## 
## Moran I statistic standard deviate = 10.391, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.537741479      -0.031250000       0.002998588

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=5)",
                  "Radial Distance Weight (dmax=105)",
                  "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.6031197 1.952469e-12
## 2     Radial Distance Weight (dmax=105)    0.7129655 8.132416e-22
## 3      Invers Distance Weight (alpha=1)    0.2765689 2.208934e-13
## 4      Invers Distance Weight (alpha=2)    0.5846982 2.573622e-09
## 5 Exponential Distance Weight (alpha=1)    0.2656131 4.911509e-27
## 6 Exponential Distance Weight (alpha=2)    0.5377415 1.365822e-25
## 7                       Rook Contiguity    0.7677398 2.604196e-10
## 8                      Queen Contiguity    0.7677398 2.604196e-10

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

##                           Matriks Bobot Indeks Moran      p-value
## 1              K-Nearest Neighbor (k=5)    0.6031197 1.952469e-12
## 2     Radial Distance Weight (dmax=105)    0.7129655 8.132416e-22
## 3      Invers Distance Weight (alpha=1)    0.2765689 2.208934e-13
## 4      Invers Distance Weight (alpha=2)    0.5846982 2.573622e-09
## 5 Exponential Distance Weight (alpha=1)    0.2656131 4.911509e-27
## 6 Exponential Distance Weight (alpha=2)    0.5377415 1.365822e-25
## 7                       Rook Contiguity    0.7677398 2.604196e-10
## 8                      Queen Contiguity    0.7677398 2.604196e-10

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

Matriks pembobot spasial yang memiliki Indeks Moran tertinggi adalah rook contiguity dan queen contiguity 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 queen contiguity.

Metode OLS

ols <- lm(PPM ~ RLS + PPK + UHH + HLS, data = data)
summary(ols)

## 
## Call:
## lm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3129 -1.9354 -0.2871  0.9077  5.9370 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -6.847e+01  4.381e+01  -1.563  0.12928   
## RLS         -2.401e+00  7.186e-01  -3.341  0.00238 **
## PPK         -7.114e-04  3.568e-04  -1.994  0.05597 . 
## UHH          1.241e+00  5.764e-01   2.153  0.04011 * 
## HLS          1.321e+00  1.146e+00   1.153  0.25851   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.578 on 28 degrees of freedom
## Multiple R-squared:  0.6691, Adjusted R-squared:  0.6219 
## F-statistic: 14.16 on 4 and 28 DF,  p-value: 1.954e-06

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

## [1] "Nilai AIC OLS 162.724559457354"

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.80478, p-value = 0.0333

Karena p-value < 5%, maka tolak H0. Artinya sisaan tidak 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 = 3.3554, df = 4, p-value = 0.5002

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 > 5.

car::vif(ols)

##      RLS      PPK      UHH      HLS 
## 4.478456 2.704217 1.834623 1.667146

Berdasarkan perhitungan, diperoleh nilai VIF < 5 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.6761, p-value = 0.1009
## alternative hypothesis: true autocorrelation is greater than 0

Karena p-value > 0.05, maka gagal tolak H0. Artinya, antargalat tidak 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.queen
lm.morantest(ols, listw=ww, alternative="two.sided")

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 0.71837, p-value = 0.4725
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.001584773     -0.088268932      0.015645027

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 = 0.24288, p-value = 0.8081
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.001584773      -0.031250000       0.018276611

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 = 6.2127, p-value = 5.208e-10
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.76773982       -0.03125000        0.01653951

moran.plot(data$PPM, ww, labels=data$`NAMA KABUPATEN`)

Uji Indeks Moran pada peubah X1/RLS

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

## 
##  Moran I test under randomisation
## 
## data:  data$RLS  
## weights: ww    
## 
## Moran I statistic standard deviate = 4.2142, p-value = 2.507e-05
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.53376200       -0.03125000        0.01797592

Uji Indeks Moran pada peubah X2/PPK

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

## 
##  Moran I test under randomisation
## 
## data:  data$PPK  
## weights: ww    
## 
## Moran I statistic standard deviate = 5.4395, p-value = 5.342e-08
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.70263327       -0.03125000        0.01820257

Uji Indeks Moran pada peubah X3/UHH

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

## 
##  Moran I test under randomisation
## 
## data:  data$UHH  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.5232, p-value = 0.0004263
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.45449156       -0.03125000        0.01900772

Uji Indeks Moran pada peubah X4/HLS

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

## 
##  Moran I test under randomisation
## 
## data:  data$HLS  
## weights: ww    
## 
## Moran I statistic standard deviate = 1.2661, p-value = 0.2055
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.1375906        -0.0312500         0.0177847

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

Peubah = c("Y", "X1", "X2","X3","X4", "Sisaan")
Indeks_Moran = c(my$estimate[1], mx1$estimate[1], mx2$estimate[1],mx3$estimate[1],mx4$estimate[1], merr$estimate[1])
p_value = c(my$p.value, mx1$p.value, mx2$p.value,mx3$p.value,mx4$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.767739817 5.208393e-10
## 2     X1  0.533762002 2.506973e-05
## 3     X2  0.702633271 5.342396e-08
## 4     X3  0.454491559 4.263306e-04
## 5     X4  0.137590564 2.054923e-01
## 6 Sisaan  0.001584773 8.081008e-01

Berdasarkan hasil pengujian Indeks Moran di atas menunjukkan bahwa nilai indeks Moran pada peubah Y/PPM, X1/RLS, X2/PPK, dan X3/UHH 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/RLS, X2/PPK, dan X3/UHH regresi klasik. Terlihat pula nilai indeks moran positif pada peubah Y/PPM, X1/RLS, X2/PPK, dan X3/UHH regresi klasik yang menunjukkan adanya autokorelasi positif.

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

Oleh karenanya, 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 ~ RLS + PPK + UHH + HLS, data = data)
## test weights: listw
##  
##           statistic parameter  p.value   
## RSerr    0.00012197         1 0.991188   
## adjRSerr 4.48276650         1 0.034238 * 
## RSlag    3.28066772         1 0.070100 . 
## adjRSlag 7.76331225         1 0.005332 **
## SARMA    7.76343422         2 0.020615 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Output menunjukkan bahwa hasil uji robust LM SAR, robust 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 = 3.3554, df = 4, p-value = 0.5002

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 ~ RLS + PPK + UHH + HLS, data = data, listw = ww, zero.policy=TRUE)
summary(sar, Nagelkerke = T)

## 
## Call:lagsarlm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data, 
##     listw = ww, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.68182 -1.52917 -0.59518  1.10309  5.61842 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -5.4235e+01  3.8386e+01 -1.4129 0.157691
## RLS         -1.8142e+00  6.5445e-01 -2.7722 0.005568
## PPK         -3.8415e-04  3.2358e-04 -1.1872 0.235159
## UHH          9.2897e-01  5.1623e-01  1.7995 0.071935
## HLS          1.0458e+00  9.8961e-01  1.0568 0.290605
## 
## Rho: 0.33589, LR test value: 4.0616, p-value: 0.043868
## Asymptotic standard error: 0.14591
##     z-value: 2.3021, p-value: 0.021328
## Wald statistic: 5.2998, p-value: 0.021328
## 
## Log likelihood: -73.33146 for lag model
## ML residual variance (sigma squared): 4.8257, (sigma: 2.1968)
## Nagelkerke pseudo-R-squared: 0.70744 
## Number of observations: 33 
## Number of parameters estimated: 7 
## AIC: 160.66, (AIC for lm: 162.72)
## LM test for residual autocorrelation
## test value: 2.5377, p-value: 0.11115

Output di atas memperlihatkan bahwa koefisien Rho pada model SAR signifikan, dengan nilai AIC sebesar 160.66. Selain itu, terlihat pula hasil uji autokorelasi pada sisaan model memperlihatkan nilai p-value sebesar 0.11115, 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.84895, p-value = 0.02576

Karena p-value < 5%, maka tolak H0. Artinya sisaan tidak menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 4.8419, df = 4, p-value = 0.3039

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 = -0.80849, p-value = 0.4188
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.14032929       -0.03125000        0.01820247

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

Model SEM

sem <- errorsarlm(PPM ~ RLS + PPK + UHH + HLS,data=data,listw=ww)
summary(sem)

## 
## Call:errorsarlm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data, 
##     listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.32169 -1.93480 -0.31854  0.89650  5.95662 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) -6.8711e+01  4.0519e+01 -1.6958 0.0899308
## RLS         -2.3862e+00  6.6322e-01 -3.5978 0.0003209
## PPK         -7.1239e-04  3.2944e-04 -2.1625 0.0305832
## UHH          1.2423e+00  5.3406e-01  2.3262 0.0200105
## HLS          1.3218e+00  1.0569e+00  1.2506 0.2110926
## 
## Lambda: 0.012583, LR test value: 0.00065679, p-value: 0.97955
## Asymptotic standard error: 0.21061
##     z-value: 0.059747, p-value: 0.95236
## Wald statistic: 0.0035697, p-value: 0.95236
## 
## Log likelihood: -75.36195 for error model
## ML residual variance (sigma squared): 5.6377, (sigma: 2.3744)
## Number of observations: 33 
## Number of parameters estimated: 7 
## AIC: 164.72, (AIC for lm: 162.72)

Output di atas menunjukkan bahwa koefisien Lambda tidak signifikan pada taraf nyata 5%. AIC model SEM adalah 164.72. 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.80351, p-value = 0.03355

Karena p-value < 0.05 maka tolak H0, artinya sisaan tidak menyebar normal.

Asumsi Kehomogenan Ragam

H0 : Ragam Sisaan Homogen

H1 : Ragam Sisaan Tidak Homogen

bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 3.3846, df = 4, p-value = 0.4956

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.23634, p-value = 0.8132
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0006963127     -0.0312500000      0.0182708668

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

Model SARMA/GSM

SARMA <- sacsarlm(PPM ~ RLS + PPK + UHH + HLS,data=data,ww)
summary(SARMA)

## 
## Call:sacsarlm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data, 
##     listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.25204 -1.35166 -0.44601  0.79217  4.40486 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -1.5135e+01  2.8840e+01 -0.5248  0.59972
## RLS         -1.3414e+00  6.7407e-01 -1.9900  0.04659
## PPK         -2.6345e-04  2.3587e-04 -1.1169  0.26403
## UHH          4.6192e-01  3.7003e-01  1.2483  0.21191
## HLS          9.5109e-02  8.3352e-01  0.1141  0.90915
## 
## Rho: 0.53457
## Asymptotic standard error: 0.15805
##     z-value: 3.3823, p-value: 0.00071879
## Lambda: -0.69835
## Asymptotic standard error: 0.24993
##     z-value: -2.7942, p-value: 0.0052028
## 
## LR test value: 8.0074, p-value: 0.018248
## 
## Log likelihood: -71.35858 for sac model
## ML residual variance (sigma squared): 3.5534, (sigma: 1.8851)
## Number of observations: 33 
## Number of parameters estimated: 8 
## AIC: 158.72, (AIC for lm: 162.72)

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

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.62964, p-value = 0.09234

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 = 8.1036, df = 4, p-value = 0.08786

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.52521, p-value = 0.5994
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.10224464       -0.03125000        0.01827168

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%.

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"),
           "AIC" = c(AIC(ols),AIC(sar),AIC(sem), AIC(SARMA)),"p-value dari Rho"=c("NA","0.043868","NA","0.00071879"), "p-value dari Lambda"=c("NA","NA","0.012583","0.0052028"), "Kenormalan (p-value)"=c("0.0333","0.02576","0.03355","0.09234"), "Homoskedastisitas (p-value)"=c("0.5002","0.3039","0.4956","0.08786"), "Kebebasan Sisaan (p-value)"=c("0.1009","0.4188","0.8132","0.5994"))

df

##                  Model      AIC p.value.dari.Rho p.value.dari.Lambda
## 1 OLS (Regresi KlasiK) 162.7246               NA                  NA
## 2                  SAR 160.6629         0.043868                  NA
## 3                  SEM 164.7239               NA            0.012583
## 4                SARMA 158.7172       0.00071879           0.0052028
##   Kenormalan..p.value. Homoskedastisitas..p.value. Kebebasan.Sisaan..p.value.
## 1               0.0333                      0.5002                     0.1009
## 2              0.02576                      0.3039                     0.4188
## 3              0.03355                      0.4956                     0.8132
## 4              0.09234                     0.08786                     0.5994

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

Interpretasi dan Kesimpulan

Karena Model SARMA memiliki spillover, maka perlu diinterpretasi dengan memperhatikan spillovernya. Interpretasi koefisien pada model regresi spasial dijelaskan dengan efek langsung, tidak langsung, dan total dari setiap peubah.

sum <- summary(SARMA)
sum

## 
## Call:sacsarlm(formula = PPM ~ RLS + PPK + UHH + HLS, data = data, 
##     listw = ww)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.25204 -1.35166 -0.44601  0.79217  4.40486 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) -1.5135e+01  2.8840e+01 -0.5248  0.59972
## RLS         -1.3414e+00  6.7407e-01 -1.9900  0.04659
## PPK         -2.6345e-04  2.3587e-04 -1.1169  0.26403
## UHH          4.6192e-01  3.7003e-01  1.2483  0.21191
## HLS          9.5109e-02  8.3352e-01  0.1141  0.90915
## 
## Rho: 0.53457
## Asymptotic standard error: 0.15805
##     z-value: 3.3823, p-value: 0.00071879
## Lambda: -0.69835
## Asymptotic standard error: 0.24993
##     z-value: -2.7942, p-value: 0.0052028
## 
## LR test value: 8.0074, p-value: 0.018248
## 
## Log likelihood: -71.35858 for sac model
## ML residual variance (sigma squared): 3.5534, (sigma: 1.8851)
## Number of observations: 33 
## Number of parameters estimated: 8 
## AIC: 158.72, (AIC for lm: 162.72)

# Spill over
Im = impacts(SARMA, listw = ww)
Im

## Impact measures (sac, exact):
##                  Direct      Indirect         Total
## RLS dy/dx -1.4818809449 -1.4002183899 -2.8820993348
## PPK dy/dx -0.0002910393 -0.0002750009 -0.0005660402
## UHH dy/dx  0.5102977508  0.4821765861  0.9924743369
## HLS dy/dx  0.1050694392  0.0992793392  0.2043487784

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

##         Koefisien  EfekLangsung    UmpanBalik
## RLS -1.3414063361 -1.4818809449 -1.404746e-01
## PPK -0.0002634503 -0.0002910393 -2.758901e-05
## UHH  0.4619241772  0.5102977508  4.837357e-02
## HLS  0.0951094026  0.1050694392  9.960037e-03

Koefisien = C Efek Langsung = DL Efek Umpan Balik = FB Efek Tidak Langsung (ID) = DL * (FB)^0.321 Efek Total (T) = DL + ID

Berikut ini adalah interpretasinya:

Peubah Penjelas X1/Rata - rata Lama Sekolah (RLS)

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

Pengaruh tidak langsung dari peubah X1 bernilai 0.78998, 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.78998, jika peubah lainnya tetap.

Pengaruh total dari peubah X1 bernilai -0.692, artinya apabila peningkatan X1 terjadi di seluruh wilayah maka akan menurunkan Y diseluruh wilayah dengan rata-rata penurunan sebesar 0.692. Efek umpan baliknya sebesar 0.00996.

Peubah Penjelas X2/Pengeluaran per Kapita

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

Pengaruh tidak langsung dari X2 tersebut bernilai 0.0001, 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.0001, jika peubah lainnya tetap.

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

Peubah Penjelas X3/Umur Harapan Hidup (UHH)

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

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

Pengaruh total dari peubah X3 bernilai 0.7033, artinya apabila peningkatan X3 terjadi di seluruh wilayah maka akan menurunkan Y diseluruh wilayah dengan rata-rata peningkatan sebesar 0.7033. Efek umpan baliknya sebesar -0.0000276.

Peubah Penjelas X4/Harapan Lama Sekolah

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

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

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