UAS Teknik Sampling

1 Model Spasial Global

2 Package yang dibutuhkan

library(openxlsx)
library(DT)
library(spdep)
library(rgdal)
library(raster)
library(DescTools)
library(nortest)
library(car)
library(lmtest)
library(sp)
library(spatialreg)

3 Import Data

data.jabar = read.xlsx("Jabar Data (gabung).xlsx")     # Mengimport data
datatable(data.jabar, htmltools::em('Table data'),     # Menunjukan data
          extension = 'FixedColumns',
          options = list(scrollX = TRUE, fixedColums = TRUE))  

Studi Kasus: Kemiskinan di Jawa Barat

• Y : Persentase Penduduk Miskin Tahun 2016
• X : Angka Melek Huruf Tahun 2016

Jabar2<-readOGR(dsn="petaJabar2", layer="Jabar2")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\vferd\Documents\Sherly\Campus\Semester 4\Teknik Sampling\UAS\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields

4 Eksplorasi Data

4.1 Hubungan Antar Data

plot(data.jabar$EYS2016, data.jabar$p.miskin16,        # Memplot data dengan label x dan y seperti yang tertera dan desain yang diinginkan
  xlab="Angka Melek Huruf Thn.2016", 
  ylab="Persentase Penduduk Miskin Thn.2016",
  pch=20, col="orange", cex=2)
reg.klasik = lm(p.miskin16~EYS2016, data = data.jabar) # Membuat data linear
lines.lm(reg.klasik, col=2, add=T)                     # Membuat data linear dalam bentuk garis

Plot tersebut memperlihatkan adanya pola hubungan linear negatif antara angka melek huruf terhadap persentase penduduk miskin di Jawa Barat pada tahun 2016.

4.2 Plot Persentase Penduduk Miskin tahun 2016

k=16                                                   # Sebagai parameter jumlah warna yang digunakan
colfunc <- colorRampPalette(c("green", "yellow","red"))# Memilih skala warna
color <- colfunc(k)                                    # Memilih jumlah warna

Jabar2$miskin2<- data.jabar$p.miskin16
spplot(Jabar2, "miskin2", col.regions=color)           # Melakukan plot terhadap data

Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol pada data persentase kemiskinan di kabupaten/kota di Jawa Barat. Hal ini tampak dari gradasi warna yang cenderung mengumpul, seperti pada warna merah dan oranye.

4.3 Moran Test

w<-poly2nb(Jabar2)
ww<-nb2listw(w)
moran(data.jabar$p.miskin16, ww, n=length(ww$neighbours), 
      S0=Szero(ww))

## $I
## [1] 0.3932657
## 
## $K
## [1] 2.403804

moran.test(data.jabar$p.miskin16, ww,randomisation=T,  # Melakukan Moran Test
           alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  data.jabar$p.miskin16  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.0168, p-value = 0.001277
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.3932657        -0.0400000         0.0206265

4.4 Moran Plot

moran.plot(data.jabar$p.miskin16, ww, labels=data.jabar$KABKOT)

5 Model Regresi Sederhana

5.1 OLS

reg.klasik = lm(p.miskin16~EYS2016, data = data.jabar)
err.regklasik<-residuals(reg.klasik)
summary(reg.klasik)

## 
## Call:
## lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3627 -1.9604 -0.3899  1.5238  8.1125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.5732     9.0973   4.350 0.000217 ***
## EYS2016      -2.3948     0.7209  -3.322 0.002854 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.788 on 24 degrees of freedom
## Multiple R-squared:  0.315,  Adjusted R-squared:  0.2865 
## F-statistic: 11.04 on 1 and 24 DF,  p-value: 0.002854

5.2 Model Diagnostics

5.2.1 Normality Test

ad.test(err.regklasik)

## 
##  Anderson-Darling normality test
## 
## data:  err.regklasik
## A = 0.39412, p-value = 0.3495

hist(err.regklasik)

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

H0: galat model menyebar normal
H1: galat model tidak menyebar normal

5.2.2 Heteroscedastics

bptest(reg.klasik)

## 
##  studentized Breusch-Pagan test
## 
## data:  reg.klasik
## BP = 1.0644, df = 1, p-value = 0.3022

H0: ragam galat homogen
H1: ragam galat tidak homogen

5.3 Eksplorasi Auto korelasi Spatial

Uji kebebasan sisaan pada data spasial dapat dilakukan dengan uji moran menggunakan fungsi berikut.

w<-poly2nb(Jabar2)
ww<-nb2listw(w)
lm.morantest(reg.klasik, ww, alternative="two.sided")

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## weights: ww
## 
## Moran I statistic standard deviate = 3.4736, p-value = 0.0005135
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.44799554      -0.05028275       0.02057696

Selain menggunakan fungsi lm.morantest, uji moran dapat dilakukan menggunakan fungsi moran.test seperti yang dibahas pada modul pertemuan sebelumnya. Perbedaannya adalah pada fungsi pertama, input yang digunakan adalah objek lm, sedangkan pada fungsi kedua, yang digunakan sebagai input adalah data sisaan model.

moran.test(err.regklasik, ww,randomisation=F, alternative="two.sided")

## 
##  Moran I test under normality
## 
## data:  err.regklasik  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.4266, p-value = 0.0006112
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.44799554       -0.04000000        0.02028183

Terlihat pada output bahwa hasil kedua tes menunjukkan kesimpulan yang sama, yaitu tolak H0 yang menyatakan bahwa tidak terdapat autokorelasi pada sisaan model regresi klasik pada 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.

6 Lagrange Multiplier Test

LM<-lm.LMtests(reg.klasik, nb2listw(w, style="W"),
               test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
summary(LM)

##  Lagrange multiplier diagnostics for spatial dependence
## data:  
## model: lm(formula = p.miskin16 ~ EYS2016, data = data.jabar)
## weights: nb2listw(w, style = "W")
##  
##        statistic parameter  p.value   
## LMerr   8.375484         1 0.003803 **
## LMlag   6.642323         1 0.009958 **
## RLMerr  1.822499         1 0.177016   
## RLMlag  0.089338         1 0.765020   
## SARMA   8.464822         2 0.014517 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Output memperlihatkan bahwa hasil uji model SEM dan SAR sama-sama signifikan pada taraf 5%. Selanjutnya, hasil uji robust keduanya ternyata sama-sama tidak signifikan. Berdasarkan skema tersebut, kita dapat mencoba kandidat model SARMA atau GSM. Namun demikian, ada pula pendapat yang menyarankan agar kita mengambil kandidat model dengan p-value terkecil, dalam hal ini adalah model SEM ( p-value = 0.003803 ).
Untuk selanjutnya, kita akan mencoba ketiga model, SEM, SAR, dan SARMA, meskipun pada prakteknya, Anda hanya perlu memodelkan yang menurut Anda terbaik saja.

7 Model Regresi Spasial

7.1 Model SEM

sem<-errorsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(sem)

## 
## Call:errorsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, 
##     listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.24699 -1.36122 -0.13809  1.15579  7.03019 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 36.88515    7.80246  4.7274 2.274e-06
## EYS2016     -2.17498    0.60042 -3.6224 0.0002919
## 
## Lambda: 0.61793, LR test value: 8.7676, p-value: 0.0030663
## Asymptotic standard error: 0.1576
##     z-value: 3.9208, p-value: 8.8267e-05
## Wald statistic: 15.372, p-value: 8.8267e-05
## 
## Log likelihood: -58.12466 for error model
## ML residual variance (sigma squared): 4.5459, (sigma: 2.1321)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 124.25, (AIC for lm: 131.02)

Output di atas menunjukkan bahwa koefisien Lambda signifikan pada taraf nyata 5% ( p-value = 0.0030663 ). AIC model SEM adalah sebesar 124.25. Selanjutnya kita akan coba memeriksa sisaan model SEM ini.

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

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

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

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.096995, p-value = 0.9227
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.05288671       -0.04000000        0.01765155

Terlihat pada output di atas bahwa sisaan telah memenuhi asumsi kenormalan, kehomogenan ragam, serta kebebasan.

7.2 Model SAR

sar<-lagsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(sar)

## 
## Call:
## lagsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.17670 -0.93185 -0.11318  0.91353  7.69131 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 28.22898    7.66250  3.6840 0.0002296
## EYS2016     -1.94997    0.57325 -3.4016 0.0006700
## 
## Rho: 0.59078, LR test value: 7.9343, p-value: 0.0048507
## Asymptotic standard error: 0.1559
##     z-value: 3.7894, p-value: 0.00015101
## Wald statistic: 14.36, p-value: 0.00015101
## 
## Log likelihood: -58.54132 for lag model
## ML residual variance (sigma squared): 4.7513, (sigma: 2.1798)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 125.08, (AIC for lm: 131.02)
## LM test for residual autocorrelation
## test value: 0.036687, p-value: 0.8481

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

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

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

Berdasarkan output di atas, pada taraf 5% dapat disimpulkan bahwa sisaan model telah memenuhi asumsi kenormalan dan kehomogenan ragam.

7.3 Model GSM/SARMA

gsm<-sacsarlm(p.miskin16~EYS2016,data=data.jabar,nb2listw(w))
summary(gsm)

## 
## Call:
## sacsarlm(formula = p.miskin16 ~ EYS2016, data = data.jabar, listw = nb2listw(w))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0817 -1.2815 -0.1627  1.1317  6.7128 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 38.02312    8.19667  4.6388 3.504e-06
## EYS2016     -2.15311    0.59378 -3.6261 0.0002877
## 
## Rho: -0.14689
## Asymptotic standard error: 0.38983
##     z-value: -0.3768, p-value: 0.70633
## Lambda: 0.69352
## Asymptotic standard error: 0.23014
##     z-value: 3.0135, p-value: 0.0025829
## 
## LR test value: 8.8931, p-value: 0.011719
## 
## Log likelihood: -58.06189 for sac model
## ML residual variance (sigma squared): 4.3254, (sigma: 2.0797)
## Number of observations: 26 
## Number of parameters estimated: 5 
## AIC: 126.12, (AIC for lm: 131.02)

Output di atas memperlihatkan bahwa hanya salah satu koefisien dependensi spasial yang signifikan, yaitu Lambda. AIC model SARMA adalah sebesar 126.12.

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

## 
##  Anderson-Darling normality test
## 
## data:  err.gsm
## A = 0.55307, p-value = 0.1388

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

## 
##  Moran I test under randomisation
## 
## data:  err.gsm  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.0061095, p-value = 0.9951
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.04081914       -0.04000000        0.01797655

Berdasarkan output di atas, terlihat bahwa sisaan model SARMA telah memenuhi asumsi kenormalan, kehomogenan ragam, dan kebebasan.

8 Goodness Of Fit

Ilustrasi pada kasus ini memperlihatkan bahwa ternyata SEM merupakan model terbaik berdasarkan nilai AIC-nya. Hal ini ternyata sejalan dengan p-value nya yang juga terkecil pada uji LM.