Spasial autoregressive quantile regression
2023-06-01
PACKAGE
library(terra)## terra 1.8.54library(sp)
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(remotes)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:terra':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(quantreg)## Loading required package: SparseMlibrary(readxl)
library(sp)
library(spData)## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spdep)
library(car)## Loading required package: carDatalibrary(skedastic)## Registered S3 methods overwritten by 'registry':
## method from
## print.registry_field proxy
## print.registry_entry proxylibrary(spatialreg)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:SparseM':
##
## det##
## 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.ZeroPolicyOptionlibrary(nortest)
library(ggplot2)
library(corrplot)## corrplot 0.95 loadedIMPORT DATA
data <- data_jateng <- read_excel("C:/Users/halfi/Downloads/SAR SARQR BU RINA/SAR SARQR BU RINA/data.xlsx")## New names:
## • `` -> `...9`head(data)## # A tibble: 6 × 9
## KABKOT Y X1 X2 X3 X4 X5 X6 ...9
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Kabupaten Cilacap 71.0 5.37 190960 9.62 65.6 327 402074 1
## 2 Kabupaten Banyumas 73.2 4.88 220470 6.05 64.8 389 475654 2
## 3 Kabupaten Purbalingga 69.5 2.75 145330 5.23 73.5 160 413194 3
## 4 Kabupaten Banjarnegara 68.6 2.8 141250 6.38 72.2 152 432505 4
## 5 Kabupaten Kebumen 70.8 3.72 196160 5.92 71.6 164 483725 5
## 6 Kabupaten Purworejo 73.6 2.1 82640 4.45 72.4 150 449772 6str(data)## tibble [35 × 9] (S3: tbl_df/tbl/data.frame)
## $ KABKOT: chr [1:35] "Kabupaten Cilacap" "Kabupaten Banyumas" "Kabupaten Purbalingga" "Kabupaten Banjarnegara" ...
## $ Y : num [1:35] 71 73.2 69.5 68.6 70.8 ...
## $ X1 : num [1:35] 5.37 4.88 2.75 2.8 3.72 2.1 2.42 3.54 2.92 3.45 ...
## $ X2 : num [1:35] 190960 220470 145330 141250 196160 ...
## $ X3 : num [1:35] 9.62 6.05 5.23 6.38 5.92 4.45 5.01 4.97 4.92 4.31 ...
## $ X4 : num [1:35] 65.7 64.8 73.5 72.2 71.6 ...
## $ X5 : num [1:35] 327 389 160 152 164 150 118 284 299 366 ...
## $ X6 : num [1:35] 402074 475654 413194 432505 483725 ...
## $ ...9 : num [1:35] 1 2 3 4 5 6 7 8 9 10 ...datax = data[2:8]
head(datax)## # A tibble: 6 × 7
## Y X1 X2 X3 X4 X5 X6
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 71.0 5.37 190960 9.62 65.6 327 402074
## 2 73.2 4.88 220470 6.05 64.8 389 475654
## 3 69.5 2.75 145330 5.23 73.5 160 413194
## 4 68.6 2.8 141250 6.38 72.2 152 432505
## 5 70.8 3.72 196160 5.92 71.6 164 483725
## 6 73.6 2.1 82640 4.45 72.4 150 449772attach(datax)library(sf)
library(sp)
jateng_sf = st_read("C:/Users/halfi/Downloads/SAR SARQR BU RINA/SAR SARQR BU RINA/RBI_50K_2023_Jawa Tengah/RBI_50K_2023_Jawa Tengah.shp")## Reading layer `RBI_50K_2023_Jawa Tengah' from data source
## `C:\Users\halfi\Downloads\SAR SARQR BU RINA\SAR SARQR BU RINA\RBI_50K_2023_Jawa Tengah\RBI_50K_2023_Jawa Tengah.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 35 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 108.5558 ymin: -8.211806 xmax: 111.6914 ymax: -5.725371
## Geodetic CRS: WGS 84jateng = as(jateng_sf, "Spatial")
class(jateng)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"coords = coordinates(jateng)EXPLORATORY DATA
summary(datax)## Y X1 X2 X3
## Min. :67.03 Min. :0.330 Min. : 8650 Min. :1.760
## 1st Qu.:70.78 1st Qu.:2.255 1st Qu.: 75820 1st Qu.:4.355
## Median :73.15 Median :2.790 Median : 97180 Median :5.010
## Mean :73.50 Mean :2.857 Mean :109470 Mean :5.345
## 3rd Qu.:75.89 3rd Qu.:3.580 3rd Qu.:143940 3rd Qu.:6.505
## Max. :84.35 Max. :5.430 Max. :290660 Max. :9.640
## X4 X5 X6
## Min. :64.75 Min. : 89.0 Min. : 402074
## 1st Qu.:68.63 1st Qu.: 147.5 1st Qu.: 454754
## Median :70.99 Median : 179.0 Median : 522326
## Mean :71.12 Mean : 229.6 Mean : 584428
## 3rd Qu.:73.72 3rd Qu.: 237.0 3rd Qu.: 594390
## Max. :79.57 Max. :1356.0 Max. :1536037library(sp)
jateng_sp<-as(jateng_sf,"Spatial")
library(RColorBrewer)
colfunc <- colorRampPalette(c("#EEE2DC", "#DA7B93","#5D001E"))
#colorRampPalette menentukan warna apa saja yang kita inginkan
jateng_sf$Y <- Y
spplot(jateng_sp, "Y",
col.regions = colfunc(6),
cuts = 3,
at = seq(min(Y), max(Y), length.out = 4),
main = "Peta Sebaran Indeks Pembangunan Manusia di Jawa Tengah")
#KORELASI
korelasi = cor(datax)
korelasi## Y X1 X2 X3 X4 X5
## Y 1.00000000 -0.4440284 -0.63110435 -0.04075512 -0.2202671 0.44930786
## X1 -0.44402837 1.0000000 0.85455527 0.41410116 -0.2510019 0.37214764
## X2 -0.63110435 0.8545553 1.00000000 0.33424492 -0.2047442 0.05084769
## X3 -0.04075512 0.4141012 0.33424492 1.00000000 -0.4868062 0.25782272
## X4 -0.22026707 -0.2510019 -0.20474417 -0.48680618 1.0000000 -0.13671246
## X5 0.44930786 0.3721476 0.05084769 0.25782272 -0.1367125 1.00000000
## X6 0.87002747 -0.4346296 -0.56762791 0.10913585 -0.1872612 0.36429127
## X6
## Y 0.8700275
## X1 -0.4346296
## X2 -0.5676279
## X3 0.1091359
## X4 -0.1872612
## X5 0.3642913
## X6 1.0000000corrplot(korelasi,
tl.col="black",
tl.srt=360,
method = "color",
addCoef.col = "black") # Bentuk visualisasi
Analisis Regresi
reg=lm(Y~X1+X2+X3+X4+X5+X6,data=datax)
summary(reg)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = datax)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.780 -1.372 -0.201 1.022 3.123
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.216e+01 8.362e+00 11.021 1.08e-11 ***
## X1 -4.239e-01 6.424e-01 -0.660 0.51467
## X2 -1.619e-05 1.142e-05 -1.417 0.16751
## X3 -3.909e-01 1.919e-01 -2.036 0.05126 .
## X4 -2.975e-01 1.028e-01 -2.894 0.00729 **
## X5 6.655e-03 2.186e-03 3.044 0.00504 **
## X6 1.035e-05 2.028e-06 5.102 2.10e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.744 on 28 degrees of freedom
## Multiple R-squared: 0.8739, Adjusted R-squared: 0.8469
## F-statistic: 32.34 on 6 and 28 DF, p-value: 2.401e-11#Prediksi Presentase Kemiskinan
Yh_Reg=reg$coefficients[1]+(reg$coefficients[2]*X1)+(reg$coefficients[3]*X2)+(reg$coefficients[4]*X3)+(reg$coefficients[5]*X4)+(reg$coefficients[6]*X5)+(reg$coefficients[7]*X6)Uji Asumsi
#Uji Multikolinearitas
vif(reg)## X1 X2 X3 X4 X5 X6
## 7.224138 5.453935 1.592232 1.438354 2.398385 2.601139#Uji Normalitas
ressi=residuals(reg)
lillie.test(ressi)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: ressi
## D = 0.11492, p-value = 0.2858#Uji Independensi
dwtest(reg)##
## Durbin-Watson test
##
## data: reg
## DW = 1.4528, p-value = 0.0238
## alternative hypothesis: true autocorrelation is greater than 0#Uji Identik
bptest(reg)##
## studentized Breusch-Pagan test
##
## data: reg
## BP = 4.4947, df = 6, p-value = 0.61REGRESI KUANTIL MULTIVARIATE
model_kuantil=rq(Y~X1+X2+X3+X4+X5+X6,data=datax,tau=c(0.1,0.25,0.5,0.75,0.9))
summary(model_kuantil,se="iid",method="br")## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = datax)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 88.32905 9.08976 9.71743 0.00000
## X1 -0.80696 0.69822 -1.15573 0.25756
## X2 -0.00001 0.00001 -0.62179 0.53911
## X3 -0.29162 0.20863 -1.39778 0.17316
## X4 -0.28381 0.11176 -2.53938 0.01694
## X5 0.00647 0.00238 2.72313 0.01101
## X6 0.00001 0.00000 5.23435 0.00001
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = datax)
##
## tau: [1] 0.25
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 99.90583 10.15586 9.83726 0.00000
## X1 -0.64146 0.78011 -0.82227 0.41787
## X2 -0.00001 0.00001 -0.60386 0.55080
## X3 -0.65244 0.23310 -2.79898 0.00918
## X4 -0.41603 0.12487 -3.33168 0.00243
## X5 0.00651 0.00266 2.45359 0.02063
## X6 0.00001 0.00000 4.61592 0.00008
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = datax)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 103.92402 12.70484 8.17987 0.00000
## X1 -0.55548 0.97591 -0.56919 0.57377
## X2 -0.00001 0.00002 -0.78627 0.43832
## X3 -0.57634 0.29160 -1.97647 0.05803
## X4 -0.44865 0.15621 -2.87201 0.00769
## X5 0.00603 0.00332 1.81591 0.08011
## X6 0.00001 0.00000 3.25029 0.00300
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = datax)
##
## tau: [1] 0.75
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 84.54457 6.12181 13.81039 0.00000
## X1 -0.60741 0.47024 -1.29171 0.20702
## X2 -0.00003 0.00001 -3.13288 0.00403
## X3 -0.18929 0.14051 -1.34721 0.18871
## X4 -0.17381 0.07527 -2.30908 0.02853
## X5 0.01522 0.00160 9.50830 0.00000
## X6 0.00001 0.00000 5.44515 0.00001
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = datax)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 90.21812 7.93940 11.36334 0.00000
## X1 -0.48382 0.60986 -0.79333 0.43426
## X2 -0.00003 0.00001 -2.33538 0.02692
## X3 -0.52273 0.18223 -2.86856 0.00776
## X4 -0.21637 0.09762 -2.21648 0.03496
## X5 0.01365 0.00208 6.57474 0.00000
## X6 0.00001 0.00000 3.95346 0.00048#Prediksi Kemiskinan setiap Kuantil
Yhat_1=model_kuantil$coefficients[1,1]+(model_kuantil$coefficients[2,1]*data$X1)+(model_kuantil$coefficients[3,1]*data$X2)+(model_kuantil$coefficients[4,1]*data$X3)+(model_kuantil$coefficients[5,1]*data$X4)+(model_kuantil$coefficients[6,1]*data$X5)+(model_kuantil$coefficients[7,1]*data$X6)
Yhat_2=model_kuantil$coefficients[1,2]+(model_kuantil$coefficients[2,2]*data$X1)+(model_kuantil$coefficients[3,2]*data$X2)+(model_kuantil$coefficients[4,2]*data$X3)+(model_kuantil$coefficients[5,2]*data$X4)+(model_kuantil$coefficients[6,2]*data$X5)+(model_kuantil$coefficients[7,2]*data$X6)
Yhat_3=model_kuantil$coefficients[1,3]+(model_kuantil$coefficients[2,3]*data$X1)+(model_kuantil$coefficients[3,3]*data$X2)+(model_kuantil$coefficients[4,3]*data$X3)+(model_kuantil$coefficients[5,3]*data$X4)+(model_kuantil$coefficients[6,3]*data$X5)+(model_kuantil$coefficients[7,3]*data$X6)
Yhat_4=model_kuantil$coefficients[1,4]+(model_kuantil$coefficients[2,4]*data$X1)+(model_kuantil$coefficients[3,4]*data$X2)+(model_kuantil$coefficients[4,4]*data$X3)+(model_kuantil$coefficients[5,4]*data$X4)+(model_kuantil$coefficients[6,4]*data$X5)+(model_kuantil$coefficients[7,4]*data$X6)
Yhat_5=model_kuantil$coefficients[1,5]+(model_kuantil$coefficients[2,5]*data$X1)+(model_kuantil$coefficients[3,5]*data$X2)+(model_kuantil$coefficients[4,5]*data$X3)+(model_kuantil$coefficients[5,5]*data$X4)+(model_kuantil$coefficients[6,5]*data$X5)+(model_kuantil$coefficients[7,5]*data$X6)PLOT REGRESI KUANTIL UNIVARIATE
plot(X1,Y,type="n",xlab="Persentase Jumlah Penduduk", ylab="Indeks Pembangunan Manusia")
points(X1,Y,cex=1.4,col="blue")
abline(rq(Y~X1, tau=0.1), col="red")
abline(rq(Y~X1, tau=0.25), col="orange")
abline(rq(Y~X1, tau=0.5), col="black")
abline(rq(Y~X1, tau=0.75), col="green")
abline(rq(Y~X1, tau=0.9), col="violet")
legend(4.55,84.6,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
plot(X2,Y,type="n",xlab="Jumlah Penduduk Miskin", ylab="Indeks Pembangunan Manusia")
points(X2,Y,cex=1.4,col="blue")
abline(rq(Y~X2, tau=0.1), col="red")
abline(rq(Y~X2, tau=0.25), col="orange")
abline(rq(Y~X2, tau=0.5), col="black")
abline(rq(Y~X2, tau=0.75), col="green")
abline(rq(Y~X2, tau=0.9), col="violet")
legend(242200,84.6,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
plot(X3,Y,type="n",xlab="Tingkat Pengangguran Terbuka", ylab="Indeks Pembangunan Manusia")
points(X3,Y,cex=1.4,col="blue")
abline(rq(Y~X3, tau=0.1), col="red")
abline(rq(Y~X3, tau=0.25), col="orange")
abline(rq(Y~X3, tau=0.5), col="black")
abline(rq(Y~X3, tau=0.75), col="green")
abline(rq(Y~X3, tau=0.9), col="violet")
legend(8.25,84.6,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
plot(X4,Y,type="n",xlab="Tingkat Partisipasi Angkatan Kerja", ylab="Indeks Pembangunan Manusia")
points(X4,Y,cex=1.4,col="blue")
abline(rq(Y~X4, tau=0.1), col="red")
abline(rq(Y~X4, tau=0.25), col="orange")
abline(rq(Y~X4, tau=0.5), col="black")
abline(rq(Y~X4, tau=0.75), col="green")
abline(rq(Y~X4, tau=0.9), col="violet")
legend(77,84.6,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
plot(X5,Y,type="n",xlab="Kriminalitas", ylab="Indeks Pembangunan Manusia")
points(X5,Y,cex=1.4,col="blue")
abline(rq(Y~X5, tau=0.1), col="red")
abline(rq(Y~X5, tau=0.25), col="orange")
abline(rq(Y~X5, tau=0.5), col="black")
abline(rq(Y~X5, tau=0.75), col="green")
abline(rq(Y~X5, tau=0.9), col="violet")
legend(10.55,84.6,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
plot(X6,Y,type="n",xlab="Persentase Pengeluaran per Kapita Bukan Makanan", ylab="Indeks Pembangunan Manusia")
points(X6,Y,cex=1.4,col="blue")
abline(rq(Y~X6, tau=0.1), col="red")
abline(rq(Y~X6, tau=0.25), col="orange")
abline(rq(Y~X6, tau=0.5), col="black")
abline(rq(Y~X6, tau=0.75), col="green")
abline(rq(Y~X6, tau=0.9), col="violet")
legend(1340000,75.5,c("Q0.1","Q0.25","Q0.50","Q0.75","Q0.9"),lty=c(1,1),lwd=c(1),col=c("red","orange","black","green","violet"))
PEMBOBOTAN
#Queen Contiguity
queenw = poly2nb(jateng, queen = TRUE)
queenw1 = nb2listw(queenw, style = "W", zero.policy = TRUE)
queenw2 = nb2mat(queenw, style = "W", zero.policy = TRUE)
queenw1## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 35
## Number of nonzero links: 148
## Percentage nonzero weights: 12.08163
## Average number of links: 4.228571
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 35 1225 35 18.64242 151.0178yy = data.frame(queenw2)SPATIAL AUTOREGRESSIVE
#Moran Test
lm.morantest(reg, queenw1, alternative = "two.sided")##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = datax)
## weights: queenw1
##
## Moran I statistic standard deviate = -0.27443, p-value = 0.7838
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## -0.05017725 -0.01768544 0.01401784#Lagrange Multiplier Lag
lm.LMtests(reg, queenw1, test = c("LMlag"))## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = datax)
## test weights: listw
##
## RSlag = 0.37156, df = 1, p-value = 0.5422#Pemodelan SAR
SAR=lagsarlm(Y~X1+X2+X3+X4+X5+X6, listw=queenw1, data=datax, type="lag")
summary(SAR)##
## Call:lagsarlm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = datax,
## listw = queenw1, type = "lag")
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.64295 -1.42575 -0.09703 1.02463 3.10880
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.6518e+01 1.2639e+01 6.8452 7.637e-12
## X1 -3.6902e-01 5.7424e-01 -0.6426 0.5204720
## X2 -1.6395e-05 1.0190e-05 -1.6090 0.1076161
## X3 -3.8353e-01 1.7082e-01 -2.2453 0.0247507
## X4 -2.8248e-01 9.4980e-02 -2.9741 0.0029386
## X5 6.4662e-03 1.9548e-03 3.3078 0.0009402
## X6 1.0456e-05 1.8188e-06 5.7488 8.986e-09
##
## Rho: 0.059565, LR test value: 0.32684, p-value: 0.56753
## Asymptotic standard error: 0.11033
## z-value: 0.53988, p-value: 0.58928
## Wald statistic: 0.29147, p-value: 0.58928
##
## Log likelihood: -65.05588 for lag model
## ML residual variance (sigma squared): 2.408, (sigma: 1.5518)
## Number of observations: 35
## Number of parameters estimated: 9
## AIC: 148.11, (AIC for lm: 146.44)
## LM test for residual autocorrelation
## test value: 0.74872, p-value: 0.38688#Prediksi Persentase Penduduk Miskin
Yh_SAR=SAR$coefficients[1]+(0.39668*(queenw2%*%Y))+(SAR$coefficients[2]*X1)+(SAR$coefficients[3]*X2)+(SAR$coefficients[4]*X3)+(SAR$coefficients[5]*X4)+(SAR$coefficients[6]*X5)+(SAR$coefficients[7]*X6)#Peta Persebaran Prediksi Presentase Penduduk Miskin
colfunc <- colorRampPalette(c("#EEE2DC", "#DA7B93","#5D001E"))
#colorRampPalette menentukan warna apa saja yang kita inginkan
jateng$Yh_SAR <- Yh_SAR
spplot(jateng, "Yh_SAR",
col.regions = colfunc(6),
cuts = 3,
main = "Peta Sebaran Prediksi Indeks Pembangunan Manusia di Jawa Tengah")