Install semua packages yang di perlukan
library(sf)
library(spdep)
library(raster)
library(spatial)
library(nortest)
library(car)
library(DescTools)
library(lmtest)
library(spdep)
library(spatialreg)
library(leaflet)
library(arules)
library(readr)
Langkah pertama adalah membaca data numerik dari file CSV menggunakan fungsi read_csv() dari paket readr. File CSV tersebut berisi data tentang provinsi Jawa Barat.
#input data numerik
data.jabar <- read_csv("C:/Users/acer/Downloads/Jabar (1).csv")
head(data.jabar)
## # A tibble: 6 × 24
## p.miskin15 p.miskin16 j.miskin15 j.miskin16 AHH2015 AHH2016 EYS2015 EYS2016
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 8.96 8.83 487. 491. 70.6 70.6 11.8 12.0
## 2 8.96 8.13 218. 199. 70.0 70.1 12.1 12.2
## 3 12.2 11.6 274. 261. 69.3 69.4 11.8 11.9
## 4 8.00 7.61 281. 273. 73.1 73.1 12.1 12.4
## 5 12.8 11.6 326. 299. 70.7 70.8 11.6 11.7
## 6 12.0 11.2 208. 196. 68.4 68.5 12.4 12.5
## # ℹ 16 more variables: MYS2015 <dbl>, MYS2016 <dbl>, EXP2015 <dbl>,
## # EXP2016 <dbl>, APM.SD15 <dbl>, APM.SMP15 <dbl>, APM.SMA15 <dbl>,
## # APM.PT15 <dbl>, APK.SD15 <dbl>, APK.SMP15 <dbl>, APK.SMA15 <dbl>,
## # APK.PT15 <dbl>, APS.USIA15 <dbl>, APS.USIA2 <dbl>, APS.USIA3 <dbl>,
## # APS.USIA4 <dbl>
Membaca Data Spasial: Kemudian, data spasial dalam format shapefile dibaca menggunakan fungsi st_read() dari paket sf
#Input Map Data
jabar <- st_read(dsn="C:/Users/acer/Downloads",layer="Jabar2")
## Reading layer `Jabar2' from data source `C:\Users\acer\Downloads' using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS: NA
jabar1 <- jabar[jabar$PROVINSI=="JAWA BARAT",]
plot(jabar1)
Plot Persentase Penduduk Miskin: Setelah data dibagi, dilakukan plot untuk memvisualisasikan persentase penduduk miskin di Jawa Barat. Warna yang digunakan dalam plot dikonfigurasi berdasarkan jumlah grup yang telah dibagi sebelumnya.
#Devided Group
k=4
#Plot Persentase Penduduk Miskin di Jawa Barat
p.miskin.baru<-discretize(data.jabar$p.miskin15,method="interval",categories=k) #mengubah variabel kontinyu menjadi variabel kategori (faktor) yang sesuai
jabar1$miskin<-p.miskin.baru
palette(rainbow(k))
plot(jabar1, ,col=1:k)
legend("topright",levels(jabar1$miskin), pch=15,col=1:k)
Plot dan Korelasi: Dilakukan plot dari beberapa variabel numerik terhadap persentase penduduk miskin. Selain itu, juga dihitung korelasi antara variabel-variabel tersebut dengan persentase penduduk miskin.
#Plot
par(mfrow=c(1,2))
plot(data.jabar$APS.USIA3, data.jabar$p.miskin15)
plot(data.jabar$EYS2016, data.jabar$p.miskin15)
par(mfrow=c(1,1))
#Correlation
cor(data.jabar$APS.USIA3, data.jabar$p.miskin15)
## [1] -0.4937358
cor(data.jabar$EYS2016, data.jabar$p.miskin15)
## [1] -0.5759777
Analisis Regresi: Dijalankan analisis regresi linier untuk memodelkan hubungan antara variabel independen tertentu dengan persentase penduduk miskin. Ini dilakukan dengan menggunakan fungsi lm().
#Regresion Analysis
reg.klasik <- lm(p.miskin15 ~ APM.SMA15+EYS2015,data.jabar)
summary(reg.klasik)
##
## Call:
## lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data = data.jabar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4010 -2.0415 -0.2338 1.6251 8.4105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.43620 13.48511 2.331 0.0289 *
## APM.SMA15 -0.07242 0.09943 -0.728 0.4737
## EYS2015 -1.38467 1.43364 -0.966 0.3442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.217 on 23 degrees of freedom
## Multiple R-squared: 0.2573, Adjusted R-squared: 0.1928
## F-statistic: 3.985 on 2 and 23 DF, p-value: 0.03266
Evaluasi Residual: Residual dari model regresi linier dievaluasi untuk memastikan bahwa mereka memenuhi asumsi-asumsi dari model tersebut. Ini meliputi uji normalitas, uji homoskedastisitas, dan uji independensi.
#Call Residual
err.regklasik <- residuals(reg.klasik)
#Normality of Residual
ad.test(err.regklasik)
##
## Anderson-Darling normality test
##
## data: err.regklasik
## A = 0.26934, p-value = 0.6511
qqnorm(err.regklasik,datax=T)
qqline(rnorm(length(err.regklasik),mean(err.regklasik),sd(err.regklasik)),datax=T, col="blue")
#Independence of Residual
RunsTest(err.regklasik)
##
## Runs Test for Randomness
##
## data: err.regklasik
## runs = 13, m = 13, n = 13, p-value = 0.8358
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x)
## -0.2338376
#Homogeneous variance of Residual
bptest(reg.klasik)
##
## studentized Breusch-Pagan test
##
## data: reg.klasik
## BP = 1.7501, df = 2, p-value = 0.4168
Uji Moran: Uji Moran digunakan untuk mengevaluasi apakah terdapat pola spasial (clustering) dalam persentase penduduk miskin di wilayah Jawa Barat.
#Weight Matrix
w<-poly2nb(jabar) # fungsi untuk menghitung dan menampilkan data relasi dari satu kabupaten dengan wilayah kabupaten lain yang bersinggungan (merubah polygon menjadi daftar tetangga)
ww<-nb2listw(w) # merubah perintah neighbor menjadi matriks daftar
#Indeks Moran
moran(data.jabar$p.miskin15, ww, n=length(ww$neighbours), S0=Szero(ww))
## $I
## [1] 0.4218876
##
## $K
## [1] 2.259555
#Moran Test
moran.test(data.jabar$p.miskin15, ww,randomisation=T, alternative="greater")
##
## Moran I test under randomisation
##
## data: data.jabar$p.miskin15
## weights: ww
##
## Moran I statistic standard deviate = 3.2057, p-value = 0.0006736
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.42188758 -0.04000000 0.02075945
#Moran Plot
moran.plot(data.jabar$p.miskin15, ww, labels=data.jabar$KABKOT)
#Indeks Moran Residual
moran(err.regklasik, ww, n=length(ww$neighbours), S0=Szero(ww))
## $I
## [1] 0.445957
##
## $K
## [1] 3.583428
#Moran Test
moran.test(err.regklasik, ww,randomisation=T, alternative="greater")
##
## Moran I test under randomisation
##
## data: err.regklasik
## weights: ww
##
## Moran I statistic standard deviate = 3.4765, p-value = 0.000254
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.44595701 -0.04000000 0.01953931
Langrange Multiplier: Dilakukan pengujian dengan menggunakan Langrange Multiplier untuk mengetahui adanya efek spasial dalam model regresi.
#Langrange Multiplier
LM <- lm.RStests(reg.klasik, nb2listw(w, style="W"),test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
LM
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data =
## data.jabar)
## test weights: nb2listw(w, style = "W")
##
## RSerr = 8.2994, df = 1, p-value = 0.003966
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data =
## data.jabar)
## test weights: nb2listw(w, style = "W")
##
## RSlag = 7.4623, df = 1, p-value = 0.0063
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data =
## data.jabar)
## test weights: nb2listw(w, style = "W")
##
## adjRSerr = 0.83909, df = 1, p-value = 0.3597
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data =
## data.jabar)
## test weights: nb2listw(w, style = "W")
##
## adjRSlag = 0.0019844, df = 1, p-value = 0.9645
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data =
## data.jabar)
## test weights: nb2listw(w, style = "W")
##
## SARMA = 8.3014, df = 2, p-value = 0.01575
Model Spasial: Dibuat beberapa model spasial, seperti Spatial Autoregressive (SAR) Model, Spatial Error Model (SEM), dan SARMA Model. Setiap model dievaluasi dan diuji untuk mengecek apakah mereka memenuhi asumsi-asumsi yang diperlukan.
#SAR Model
sar <- lagsarlm(p.miskin15 ~ APM.SMA15+EYS2015,data = data.jabar,nb2listw(w))
summary(sar)
##
## Call:lagsarlm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data = data.jabar,
## listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.280080 -1.497115 -0.024717 1.223175 7.933720
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 25.8953761 10.2861202 2.5175 0.01182
## APM.SMA15 -0.0039211 0.0735545 -0.0533 0.95749
## EYS2015 -1.7941629 1.0620783 -1.6893 0.09116
##
## Rho: 0.64001, LR test value: 9.1582, p-value: 0.002476
## Asymptotic standard error: 0.14467
## z-value: 4.4239, p-value: 9.6929e-06
## Wald statistic: 19.571, p-value: 9.6929e-06
##
## Log likelihood: -61.10011 for lag model
## ML residual variance (sigma squared): 5.6547, (sigma: 2.378)
## Number of observations: 26
## Number of parameters estimated: 5
## AIC: 132.2, (AIC for lm: 139.36)
## LM test for residual autocorrelation
## test value: 0.00039203, p-value: 0.9842
#Call Residual
err.sar<-residuals(sar)
#Normality of Residual
ad.test(err.sar)
##
## Anderson-Darling normality test
##
## data: err.sar
## A = 0.50727, p-value = 0.1825
#Homogeneous Variance of Residual
bptest.Sarlm(sar)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 1.1644, df = 2, p-value = 0.5587
#Independence of Residual
RunsTest(err.sar)
##
## Runs Test for Randomness
##
## data: err.sar
## runs = 15, m = 13, n = 13, p-value = 0.8358
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x)
## -0.02471711
#SEM Model
sem <- errorsarlm(p.miskin15 ~ APM.SMA15+EYS2015, data = data.jabar, nb2listw(w))
summary(sem)
##
## Call:errorsarlm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data = data.jabar,
## listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.356004 -1.139456 0.047083 1.137132 6.885408
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 39.432014 9.267654 4.2548 2.092e-05
## APM.SMA15 0.045494 0.064742 0.7027 0.482249
## EYS2015 -2.600407 0.916477 -2.8374 0.004548
##
## Lambda: 0.70322, LR test value: 11.021, p-value: 0.00090102
## Asymptotic standard error: 0.13335
## z-value: 5.2735, p-value: 1.3381e-07
## Wald statistic: 27.81, p-value: 1.3381e-07
##
## Log likelihood: -60.16888 for error model
## ML residual variance (sigma squared): 5.0838, (sigma: 2.2547)
## Number of observations: 26
## Number of parameters estimated: 5
## AIC: 130.34, (AIC for lm: 139.36)
#Call Residual
err.sem<-residuals(sem)
#Normality of Residual
ad.test(err.sem)
##
## Anderson-Darling normality test
##
## data: err.sem
## A = 0.47144, p-value = 0.2251
#Homogeneous Variance of Residual
bptest.Sarlm(sem)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 0.64848, df = 2, p-value = 0.7231
#Independence of Residual
RunsTest(err.sem)
##
## Runs Test for Randomness
##
## data: err.sem
## runs = 17, m = 13, n = 13, p-value = 0.3131
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x)
## 0.04708313
#SARMA Model
sarma<-sacsarlm(p.miskin15 ~ APM.SMA15+EYS2015, data = data.jabar,nb2listw(w))
summary(sarma)
##
## Call:sacsarlm(formula = p.miskin15 ~ APM.SMA15 + EYS2015, data = data.jabar,
## listw = nb2listw(w))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2875547 -1.0473692 0.0076906 1.1375755 6.7622790
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 39.678048 9.269662 4.2804 1.865e-05
## APM.SMA15 0.042042 0.066976 0.6277 0.53018
## EYS2015 -2.556479 0.958673 -2.6667 0.00766
##
## Rho: -0.060422
## Asymptotic standard error: 0.40208
## z-value: -0.15027, p-value: 0.88055
## Lambda: 0.7262
## Asymptotic standard error: 0.2184
## z-value: 3.3251, p-value: 0.00088395
##
## LR test value: 11.045, p-value: 0.0039958
##
## Log likelihood: -60.1567 for sac model
## ML residual variance (sigma squared): 5.0013, (sigma: 2.2364)
## Number of observations: 26
## Number of parameters estimated: 6
## AIC: 132.31, (AIC for lm: 139.36)
#Call Residual
err.sarma<-residuals(sarma)
#Normality of Residual
ad.test(err.sarma)
##
## Anderson-Darling normality test
##
## data: err.sarma
## A = 0.46681, p-value = 0.2312
#Homogeneous Variance of Residual
bptest.Sarlm(sarma)
##
## studentized Breusch-Pagan test
##
## data:
## BP = 0.7068, df = 2, p-value = 0.7023
#Independence of Residual
RunsTest(err.sarma)
##
## Runs Test for Randomness
##
## data: err.sarma
## runs = 17, m = 13, n = 13, p-value = 0.3131
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
## median(x)
## 0.007690637