UAS Teknik Sampling

Running Code No 4

1 Library Packages

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

2 Import Data

2.1 Data Pendidikan dan Kemiskinan, diperoleh dari BPS

data.jabar <- read.xlsx("Jabar Data (gabung).xlsx")
datatable(data.jabar)

Studi Kasus: Kemiskinan di Jawa Barat X: Persentase Penduduk miskin tahun 2016 Y: Angka Melek Huruf tahun 2016

2.2 Peta Jawa Barat (Data polygon)

jabar2 <- readOGR(dsn="D:/Documents/Campus Documents/Semester 4/Teknik Sampling/Week 16 (UAS TekSam)/petaJabar2", layer="Jabar2")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\Campus Documents\Semester 4\Teknik Sampling\Week 16 (UAS TekSam)\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields

3 Eksplorasi Data

3.1 Hubungan yang ada pada data

plot(data.jabar$EYS2016, data.jabar$p.miskin16,
     xlab="Angka Melek Huruf Thn.2016", 
     ylab="Persentase Penduduk Miskin Thn.2016",
     pch=20,
     col="blue",
     cex=2)
reg.klasik <- lm(p.miskin16~EYS2016, data = data.jabar)

lines.lm(reg.klasik, col=2, add=T)

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

3.2 Plot persentase penduduk miskin tahun 2016

k       <- 16
colfunc <- colorRampPalette(c("green","yellow","red"))
color   <- colfunc(k)

jabar2$miskin2 <- data.jabar$p.miskin16
spplot(jabar2, "miskin2", col.regions=color)

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

3.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, 
           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

3.4 Moran Plot

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

4 Classical Regression Modeling

4.1 Metode 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

4.2 Model Diagnostics

4.2.1 Uji Normalitas

\[ H_0 &= galat model menyebar normal H_a &= galat model tidak menyebar normal \]

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

4.2.2 Heteroscedastis

\[ H_0 &= ragam galat homogen H_a &= ragam galat tidak homogen \]

bptest(reg.klasik)

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

4.3 Explore for Spatial Autocorrelation

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.

5 Uji Lagrange Multiplier (LM)

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

6 Spatial Regression Modelling

6.1 Spatial Error 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

bptest.sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.38918, df = 1, p-value = 0.5327

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.

6.2 Spatial Autoregressive 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

bptest.sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.97078, df = 1, p-value = 0.3245

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

6.3 General Spatial Model (GSM)

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

bptest.sarlm(gsm)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.31097, df = 1, p-value = 0.5771

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.

Mohon diingat bahwa pada ilustrasi yang kita lakukan saat ini, kita hanya menggunakan satu peubah bebas sehingga kita tidak perlu mengkhawatirkan masalah multikolinieritas. Pada saat Anda memiliki lebih dari satu peubah bebas, pastikan Anda juga memperhatikan multikolinieritas pada model. Pemeriksaan dapat dilakukan dengan fungsi vif() pada package car.

7 Kebaikan Model

Rangkuman hasil pemodelan yang telah dilakukan sepanjang ilustrasi pada modul ini.

Rincian	OLS	SEM	SAR	GSM
AIC	131.020	124.250	125.080	126.120
p-value of Rho	NA	NA	0.005	0.706
p-value of Lambda	NA	0.003	NA	0.003
Kenormalan(p-value)	0.349	0.090	0.056	0.139
Homoskedastisitas (p-value)	0.302	0.533	0.325	0.577
Kebebasan Sisaan(p-value)	0.001	0.923	0.848	0.995

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.