Penugasan P7 ADS

Package

library(readxl)
#PACKAGE UNTUK SPASIAL
library(tmap)
library(raster)
library(spdep)
library(spatialreg)
library(rgdal)
#PACKAGE UNTUK UJI ASUMSI
library(lmtest)
library(tseries) 
library(car)
library(sf)

Load Data

spJawa = shapefile("D:/4. ANALISIS SPASIAL/Praktikum/DATA/kabkota_jawa/kabkot_jawa.shp")

Peta

tm_shape(spJawa) + tm_polygons()

Matriks Pembobot

queen.nb=poly2nb(spJawa) #Pembobot queen
queen.jawa=nb2listw(queen.nb,style="W",zero.policy=TRUE)

#Menyimpan Matriks Pembobot
bobot.queen = listw2mat(queen.jawa) #convert listw to matrix
setwd("D:/4. ANALISIS SPASIAL/Praktikum/P7")
write.csv(bobot.queen, "Matriks Bobot Queen Jawa.csv")

#modifikasi matriks, lalu load ulang
bobot=read_excel("D:/4. ANALISIS SPASIAL/Praktikum/P7/Matriks Bobot Queen Jawa.xlsx")
head(bobot)

## # A tibble: 6 x 119
##      V1    V2    V3    V4    V5    V6    V7    V8    V9   V10   V11   V12   V13
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0     0     0     0     0     1         0     0     0     0     0     0     0
## 2 0     0     0.167 0.167 0.167 0         0     0     0     0     0     0     0
## 3 0     0.167 0     0.167 0     0.167     0     0     0     0     0     0     0
## 4 0     0.25  0.25  0     0.25  0.25      0     0     0     0     0     0     0
## 5 0     0.2   0     0.2   0     0.2       0     0     0     0     0     0     0
## 6 0.167 0     0.167 0.167 0.167 0         0     0     0     0     0     0     0
## # ... with 106 more variables: V14 <dbl>, V15 <dbl>, V16 <dbl>, V17 <dbl>,
## #   V18 <dbl>, V19 <dbl>, V20 <dbl>, V21 <dbl>, V22 <dbl>, V23 <dbl>,
## #   V24 <dbl>, V25 <dbl>, V26 <dbl>, V27 <dbl>, V28 <dbl>, V29 <dbl>,
## #   V30 <dbl>, V31 <dbl>, V32 <dbl>, V33 <dbl>, V34 <dbl>, V35 <dbl>,
## #   V36 <dbl>, V37 <dbl>, V38 <dbl>, V39 <dbl>, V40 <dbl>, V41 <dbl>,
## #   V42 <dbl>, V43 <dbl>, V44 <dbl>, V45 <dbl>, V46 <dbl>, V47 <dbl>,
## #   V48 <dbl>, V49 <dbl>, V50 <dbl>, V51 <dbl>, V52 <dbl>, V53 <dbl>, ...

bobot_modif=as.matrix(bobot)
modif_bobot=mat2listw(bobot_modif,style="W",)

Moran

#MORAN TEST: PEMBOBOT QUEEN
moran.test(spJawa$Persen2021,modif_bobot,randomisation=FALSE, zero.policy = TRUE)

## 
##  Moran I test under normality
## 
## data:  spJawa$Persen2021  
## weights: modif_bobot    
## 
## Moran I statistic standard deviate = 6.542, p-value = 3.034e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.431812671      -0.008474576       0.004529460

moran.plot(spJawa$Persen2021,modif_bobot)

• Nilai Moran’s I sebesar 0.4318 mengindikasikan adanya autokorelasi spasial yang positif. Artinya, wilayah dengan persen kemiskinan tinggi bertetangga dengan wilayah dengan persen kemiskinan tinggi juga.

• Selanjutnya dilakukan pengujian dengan hipotesis sebagai berikut:

  $ H_0: I= 0 $

  $ H_1: I

  ≠

  0 $

• Didapatkan p-value sebesar 0.000 yang lebih kecil dari alfa (5%) sehingga keputusan yang diambil adalah tolak H0.

• Oleh karena itu, dengan taraf uji 5% dapat disimpulkan bahwa berdasarkan karakteristik persen kemiskinan 2021, kabupaten kota di jawa memiliki efek spasial.

Pemodelan

RLS

#REGRESI OLS
reg1=lm(spJawa$Persen2021~spJawa$RLS_2021)
summary(reg1)

## 
## Call:
## lm(formula = spJawa$Persen2021 ~ spJawa$RLS_2021)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6312 -1.9721 -0.1753  1.6121  9.1585 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      26.5430     1.4210   18.68   <2e-16 ***
## spJawa$RLS_2021  -1.9203     0.1655  -11.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.878 on 117 degrees of freedom
## Multiple R-squared:  0.5352, Adjusted R-squared:  0.5312 
## F-statistic: 134.7 on 1 and 117 DF,  p-value: < 2.2e-16

AIC_reg1 = AIC(reg1) #melihat nilai AIC
AIC_reg1

## [1] 593.2434

• Pada pemodelan menggunakan regresi linier sederhana (metode OLS), didapatkan bahwa variabel Rata Lama Sekolah signifikan berpengaruh terhadap nilai persen kemiskinan di kabupaten kota di Pulau Jawa.

#UJI ASUMSI MODEL
#NORMALITAS
res <- reg1$residuals
jarque.bera.test(res)

## 
##  Jarque Bera Test
## 
## data:  res
## X-squared = 6.4594, df = 2, p-value = 0.03957

#HOMOSKEDASTISITAS
bptest(reg1)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg1
## BP = 6.4184, df = 1, p-value = 0.01129

• Output di atas menunjukkan bahwa asumsi normalitas dan homoskedastis eror tidak terpenuhi pada tingkat signifikansi 5%

Spasial

uji_lm=lm.LMtests(reg1,modif_bobot,test = "all",zero.policy = T)
summary(uji_lm)

##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = spJawa$Persen2021 ~ spJawa$RLS_2021)
## weights: modif_bobot
##  
##        statistic parameter   p.value    
## LMerr   28.69053         1 8.492e-08 ***
## LMlag   13.50728         1 0.0002376 ***
## RLMerr  16.02325         1 6.257e-05 ***
## RLMlag   0.83999         1 0.3593990    
## SARMA   29.53052         2 3.868e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Uji LM di atas menunjukkan bahwa terdapat efek spasial

• Pada LM error, baik LM maupun Robust LM menunjukkan p-value yang lebih kecil dari 5%. Sehingga model spasial lag error (SEM) dapat digunakan

• Selain itu, pada LM lag juga menunjukkan signifikansi walaupun pada robust LM tidak signifikan. Model spasial lag dependen (SAR) mungkin dapat dipertimbangkan

• Model SARMA juga memiliki p-value yang kurang dari 5%, SARMA direkomendasikan untuk dimodelkan

SAR

sar <-lagsarlm(spJawa$Persen2021~spJawa$RLS_2021,spJawa,modif_bobot)
summary(sar,Nagelkerke=TRUE)

## 
## Call:lagsarlm(formula = spJawa$Persen2021 ~ spJawa$RLS_2021, data = spJawa, 
##     listw = modif_bobot)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.40516 -1.73925 -0.56615  1.40665  8.62973 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     20.33574    2.12648  9.5631 < 2.2e-16
## spJawa$RLS_2021 -1.59268    0.17736 -8.9801 < 2.2e-16
## 
## Rho: 0.32181, LR test value: 12.991, p-value: 0.00031305
## Asymptotic standard error: 0.084383
##     z-value: 3.8137, p-value: 0.0001369
## Wald statistic: 14.544, p-value: 0.0001369
## 
## Log likelihood: -287.1264 for lag model
## ML residual variance (sigma squared): 7.1187, (sigma: 2.6681)
## Nagelkerke pseudo-R-squared: 0.58324 
## Number of observations: 119 
## Number of parameters estimated: 4 
## AIC: 582.25, (AIC for lm: 593.24)
## LM test for residual autocorrelation
## test value: 6.8186, p-value: 0.0090214

• Model yang terbentuk yaitu\[ Persen_i=20.33-1.59RLS_i+0.32\sum w_ij*Persen_j \]

• Nilai RLS signifikan berpengaruh terhadap nilai persen kemiskinan di kabupaten/kota Pulau Jawa. Jika RLS naik 1 tahun, maka persentase kemiskinan akan turun sebesar 1,59%

• Nilai rho sebesar 0.321 dengan p-value yang kurang dari 5% menunjukkan bahwa persentase kemiskinan suatu wilayah dipengaruhi oleh persentase kemiskinan wilayah tetangganya.

• Misalnya Kota Tasikmalaya memiliki 2 tetangga, yaitu Kabupaten Ciamis dan Kabupaten Tasikmalaya. Artinya, kenaikan 1% tingkat kemiskinan Kab. Tasik akan menaikkan 0.16% tingkat kemiskinan Kota Tasik. Begitu juga jika persentase kemiskinan Kab. Ciamis naik 1%, maka persentase kemiskinan Kota Tasik akan naik 0.16%

• Nilai pseudo Rsqr yang dihasilkan sebesar 0.58. Artinya, sebesar 58% variasi persentase kemiskinan suatu wilayah dapat dijelaskan oleh variabel yang termuat dalam model.

#UJI HOMOSKEDASTISITAS
bptest.Sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 5.6868, df = 1, p-value = 0.01709

• Dalam pengujian asumsi homoskedastis, nilai p-value uji BP menunjukkan bahwa 0.01 lebih kecil dari 0.05 sehingga dapat dikatakan bahwa terdapat heteroskedastisitas eror.

SEM

sem <-errorsarlm(spJawa$Persen2021~spJawa$RLS_2021,spJawa,modif_bobot)
summary(sem,Nagelkerke=TRUE)

## 
## Call:errorsarlm(formula = spJawa$Persen2021 ~ spJawa$RLS_2021, data = spJawa, 
##     listw = modif_bobot)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -5.30213 -1.54685 -0.10254  1.50422  7.71125 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     26.37162    1.59764  16.507 < 2.2e-16
## spJawa$RLS_2021 -1.87215    0.17554 -10.665 < 2.2e-16
## 
## Lambda: 0.48234, LR test value: 23.734, p-value: 1.1062e-06
## Asymptotic standard error: 0.094264
##     z-value: 5.1169, p-value: 3.1058e-07
## Wald statistic: 26.183, p-value: 3.1058e-07
## 
## Log likelihood: -281.7548 for error model
## ML residual variance (sigma squared): 6.2805, (sigma: 2.5061)
## Nagelkerke pseudo-R-squared: 0.61922 
## Number of observations: 119 
## Number of parameters estimated: 4 
## AIC: 571.51, (AIC for lm: 593.24)

• Sehingga model yang terbentuk adalah

  \[ Persen_i=26.37-1.87RLS_i+0.48\sum w_ij*u_j \]

• Nilai RLS signifikan berpengaruh terhadap nilai persen kemiskinan di kabupaten/kota Pulau Jawa. Jika RLS naik 1 tahun, maka persentase kemiskinan akan turun sebesar 1,87%

• Nilai lambda sebesar 0.4823 dengan p-value yang kurang dari 5% menunjukkan bahwa persentase kemiskinan suatu wilayah dipengaruhi oleh eror (komponen yang tidak termuat dalam model) wilayah tetangganya.

  \[ u_T=0.48\sum 0.5u_j \]

  \[ u_T=0.48(0.5u_c+0.5u_t) \]

  \[ u_T=0.24u_c +0.24u_t \]

• Misalnya Kota Tasikmalaya (T) memiliki 2 tetangga, yaitu Kabupaten Ciamis (c) dan Kabupaten Tasikmalaya (t). Artinya, kenaikan 1 satuan komponen di luar model Kab. Tasik akan menaikkan 0.24% tingkat kemiskinan Kota Tasik. Begitu juga jika eror (komponen di luar model) Kab. Ciamis naik 1 satuan, maka persentase kemiskinan Kota Tasik akan naik 0.24%

• Nilai pseudo Rsqr yang dihasilkan sebesar 0.61. Artinya, sebesar 61% variasi persentase kemiskinan suatu wilayah dapat dijelaskan oleh variabel yang termuat dalam model.

#UJI HOMOSKEDASTISITAS
bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 7.0159, df = 1, p-value = 0.008079

• Dalam pengujian asumsi homoskedastis, nilai p-value uji BP menunjukkan bahwa 0.01 lebih kecil dari 0.05 sehingga dapat dikatakan bahwa terdapat heteroskedastisitas eror.

SARMA

sarma <-sacsarlm(spJawa$Persen2021~spJawa$RLS_2021,spJawa,modif_bobot)
summary(sarma, Nagelkerke=TRUE)

## 
## Call:sacsarlm(formula = spJawa$Persen2021 ~ spJawa$RLS_2021, data = spJawa, 
##     listw = modif_bobot)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9600 -1.5152 -0.2903  1.6004  6.5108 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     31.03104    2.26293  13.713 < 2.2e-16
## spJawa$RLS_2021 -1.84956    0.17507 -10.565 < 2.2e-16
## 
## Rho: -0.41166
## Asymptotic standard error: 0.14804
##     z-value: -2.7806, p-value: 0.0054257
## Lambda: 0.73398
## Asymptotic standard error: 0.088819
##     z-value: 8.2638, p-value: 2.2204e-16
## 
## LR test value: 28.169, p-value: 7.6422e-07
## 
## Log likelihood: -279.5373 for sac model
## ML residual variance (sigma squared): 5.2566, (sigma: 2.2927)
## Nagelkerke pseudo-R-squared: 0.63315 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 569.07, (AIC for lm: 593.24)

• Sehingga persamaan yang terbentuk adalah

  \[ Persen_i=31.03-1.84RLS_i-0.41\sum w_ij*Persen_j+0.73\sum w_ij*u_j \]

• Nilai RLS signifikan berpengaruh terhadap nilai persen kemiskinan di kabupaten/kota Pulau Jawa. Jika RLS naik 1 tahun, maka persentase kemiskinan akan turun sebesar 1,84%

• Nilai rho sebesar -0.41 dengan p-value yang kurang dari 5% menunjukkan bahwa persentase kemiskinan suatu wilayah dipengaruhi oleh persentase kemiskinan wilayah tetangganya.

• Nilai lambda sebesar 0.73 dengan p-value yang kurang dari 5% menunjukkan bahwa persentase kemiskinan suatu wilayah dipengaruhi oleh eror (komponen yang tidak termuat dalam model) wilayah tetangganya.

  \[ u_j=0.73\sum 0.5u_j \]

  \[ u_T=0.73(0.5u_c+0.5u_t) \]

  \[ u_T=0.365u_c +0.365u_t \]

• Misalnya Kota Tasikmalaya (T) memiliki 2 tetangga, yaitu Kabupaten Ciamis (c) dan Kabupaten Tasikmalaya (t). Artinya, kenaikan 1 satuan komponen di luar model Kab. Tasik akan menaikkan 0.365% tingkat kemiskinan Kota Tasik. Begitu juga jika eror (komponen di luar model) Kab. Ciamis naik 1 satuan, maka persentase kemiskinan Kota Tasik akan naik 0.365%

• Jika persentase kemiskinan Kab Tasik naik 1% maka kemiskinan Kota Tasik akan turun 0.201%. Begitu juga jika kemiskinan Kab Ciamis naik 1%, kemiskinan Kota Tasik akan turun 0.201%

• Nilai pseudo Rsqr yang dihasilkan sebesar 0.63. Artinya, sebesar 63% variasi persentase kemiskinan suatu wilayah dapat dijelaskan oleh variabel yang termuat dalam model.

#UJI HOMOSKEDASTISITAS
bptest.Sarlm(sarma)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 5.6883, df = 1, p-value = 0.01708

• Dalam pengujian asumsi homoskedastis, nilai p-value uji BP menunjukkan bahwa 0.01 lebih kecil dari 0.05 sehingga dapat dikatakan bahwa terdapat heteroskedastisitas eror.

Model Terbaik

LR.Sarlm(sem, sarma)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = -4.4351, df = 1, p-value = 0.03521
## sample estimates:
##   Log likelihood of sem Log likelihood of sarma 
##               -281.7548               -279.5373

• Pengujian di atas dilakukan untuk melihat apakah terdapat perbedaan log-likelihood model sem dan sarma.

  H0: Log likelihood kedua model sama

  H1: Terdapat perbedaan nilai log-likelihood pada kedua model

• Diperoleh p-value yang kurang dari 5% sehingga dapat dikatakan bahwa nilai log-likelihood kedua model berbeda secara signifikan.

• Log likelihood model sarma memiliki nilai yang lebih kecil dibandingkan model SEM

LR.Sarlm(sar, sarma)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = -15.178, df = 1, p-value = 9.783e-05
## sample estimates:
##   Log likelihood of sar Log likelihood of sarma 
##               -287.1264               -279.5373

• Pengujian di atas dilakukan untuk melihat apakah terdapat perbedaan log-likelihood model sar dan sarma.

  H0: Log likelihood kedua model sama

  H1: Terdapat perbedaan nilai log-likelihood pada kedua model

• Diperoleh p-value yang kurang dari 5% sehingga dapat dikatakan bahwa nilai log-likelihood kedua model berbeda secara signifikan.

• Log likelihood model sarma memiliki nilai yang lebih kecil dibandingkan model sar

  Tabel perbandingan model

  Model	AIC	Rsqr	Log Likelihood
  SAR	582.25	0.583	-287.126
  SEM	571.51	0.619	-281.755
  SARMA	569.07	0.633	-279.537

Oleh karena itu, model yang dapat digunakan adalah model SARMA