Pemodelan GWR Tingkat Kemiskinan di Indonesia
Muhammad Andriyas Musa
2024-11-07
1. Pendahuluan
Data dari Badan Pusat Statistik (BPS) menunjukkan bahwa angka kemiskinan meningkat tajam, mencapai 27,54 juta jiwa pada tahun 2021, sebelum akhirnya mengalami penurunan hingga 25,9 juta jiwa pada tahun 2023. Salah satu penyebab utama peningkatan kemiskinan ini adalah lonjakan angka pengangguran yang terjadi selama pandemi, yang meningkat dari 4,9 persen menjadi 7 persen. Selain pengangguran, Indeks Pembangunan Manusia (IPM), yang mencakup aspek kesehatan, pendidikan, dan daya beli, juga memainkan peran penting dalam menentukan kualitas hidup masyarakat. Dengan menganalisis hubungan antara tingkat pengangguran, IPM, dan angka kemiskinan, penelitian ini bertujuan untuk memberikan wawasan yang lebih dalam mengenai tantangan yang dihadapi masyarakat Indonesia dan membantu merumuskan rekomendasi kebijakan yang lebih efektif untuk mengatasi kemiskinan di masa depan.
2. Persiapan (Set-Up)
2.1. Memuat Library
library(openxlsx)
library(dplyr)
library(tidyverse)
library(gstat)
library(spgwr)
library(spdep)
library(car)
library(lmtest)
library(ggplot2)
library(sp)
library(sf)2.2. Memuat Direktori Kerja
setwd("C:/Users/Lenovo/Downloads")3. Dataset
Peta Indonesia dengan 34 Provinsi, karena data yang dipakai adalah data tahun 2022, dimana Indonesia Masih 34 Provinsi.
3.1. Mengambil Data Spasial
Untuk Data Spasial dapat diunduh di laman web GADM, dengan GADM versi 4.1 dengan 34 Provinsi.
Indonesia <- st_read("gadm41_IDN_1.json")
unique(Indonesia$NAME_1)
plot(st_geometry(Indonesia), main = "Peta Indonesia")
3.2. Mengambil Data Kasus
data = read.xlsx("dataku.xlsx")[, c("Provinsi", "PPM22", "IPM22", "TPT22")]
data[, c("Provinsi", "PPM22")] # Data Tingkat Kemiskinan## Provinsi PPM22
## 1 Aceh 14.64
## 2 SumateraUtara 8.42
## 3 SumateraBarat 5.92
## 4 Riau 6.78
## 5 Jambi 7.62
## 6 SumateraSelatan 11.90
## 7 Bengkulu 14.62
## 8 Lampung 11.57
## 9 BangkaBelitung 4.45
## 10 KepulauanRiau 6.24
## 11 JakartaRaya 4.69
## 12 JawaBarat 8.06
## 13 JawaTengah 10.93
## 14 Yogyakarta 11.34
## 15 JawaTimur 10.38
## 16 Banten 6.16
## 17 Bali 4.57
## 18 NusaTenggaraBarat 13.68
## 19 NusaTenggaraTimur 20.05
## 20 KalimantanBarat 6.73
## 21 KalimantanTengah 5.28
## 22 KalimantanSelatan 4.49
## 23 KalimantanTimur 6.31
## 24 KalimantanUtara 6.77
## 25 SulawesiUtara 7.28
## 26 SulawesiTengah 12.33
## 27 SulawesiSelatan 8.63
## 28 SulawesiTenggara 11.17
## 29 Gorontalo 15.42
## 30 SulawesiBarat 11.75
## 31 Maluku 15.97
## 32 MalukuUtara 6.23
## 33 PapuaBarat 21.33
## 34 Papua 26.56data[, c("Provinsi", "IPM22")] # Data Indeks Pembangunan Manusia## Provinsi IPM22
## 1 Aceh 72.80
## 2 SumateraUtara 72.71
## 3 SumateraBarat 73.26
## 4 Riau 73.52
## 5 Jambi 72.14
## 6 SumateraSelatan 70.90
## 7 Bengkulu 72.16
## 8 Lampung 70.45
## 9 BangkaBelitung 72.24
## 10 KepulauanRiau 76.46
## 11 JakartaRaya 81.65
## 12 JawaBarat 73.12
## 13 JawaTengah 72.79
## 14 Yogyakarta 80.64
## 15 JawaTimur 72.75
## 16 Banten 73.32
## 17 Bali 76.44
## 18 NusaTenggaraBarat 69.46
## 19 NusaTenggaraTimur 65.90
## 20 KalimantanBarat 68.63
## 21 KalimantanTengah 71.63
## 22 KalimantanSelatan 71.84
## 23 KalimantanTimur 77.44
## 24 KalimantanUtara 71.83
## 25 SulawesiUtara 73.81
## 26 SulawesiTengah 70.28
## 27 SulawesiSelatan 72.82
## 28 SulawesiTenggara 72.23
## 29 Gorontalo 69.81
## 30 SulawesiBarat 66.92
## 31 Maluku 70.22
## 32 MalukuUtara 69.47
## 33 PapuaBarat 65.89
## 34 Papua 61.39data[, c("Provinsi", "TPT22")] # Data Tingkat Pengangguran Terbuka## Provinsi TPT22
## 1 Aceh 5.97
## 2 SumateraUtara 5.47
## 3 SumateraBarat 6.17
## 4 Riau 4.40
## 5 Jambi 4.70
## 6 SumateraSelatan 4.74
## 7 Bengkulu 3.39
## 8 Lampung 4.31
## 9 BangkaBelitung 4.18
## 10 KepulauanRiau 8.02
## 11 JakartaRaya 8.00
## 12 JawaBarat 8.35
## 13 JawaTengah 5.75
## 14 Yogyakarta 3.73
## 15 JawaTimur 4.81
## 16 Banten 8.53
## 17 Bali 4.84
## 18 NusaTenggaraBarat 3.92
## 19 NusaTenggaraTimur 3.30
## 20 KalimantanBarat 4.86
## 21 KalimantanTengah 4.20
## 22 KalimantanSelatan 4.20
## 23 KalimantanTimur 6.77
## 24 KalimantanUtara 4.62
## 25 SulawesiUtara 6.51
## 26 SulawesiTengah 3.67
## 27 SulawesiSelatan 5.75
## 28 SulawesiTenggara 3.86
## 29 Gorontalo 3.25
## 30 SulawesiBarat 3.11
## 31 Maluku 6.44
## 32 MalukuUtara 4.98
## 33 PapuaBarat 5.78
## 34 Papua 3.604. Pemetaan Distribusi Data
Indonesia <- Indonesia %>%
left_join(data, by = c("NAME_1" = "Provinsi"))4.1. Pemetaan Distribusi Tingkat Kemiskinan
ggplot(Indonesia) +
geom_sf(aes(fill = PPM22), color = "white") +
scale_fill_viridis_c(option = "plasma", name = "Persentase \nPenduduk Miskin (%)") +
theme_minimal() +
labs(title = "Tingkat Kemiskinan Berdasarkan \nProvinsi Di Indonesia")
4.2. Pemetaan Distribusi Indeks Pembangunan Manusia
ggplot(Indonesia) +
geom_sf(aes(fill = IPM22), color = "white") +
scale_fill_viridis_c(option = "plasma", name = "Indeks Pembangunan\n Manusia (%)") +
theme_minimal() +
labs(title = "Indeks Pembangunan Manusia \nBerdasarkan Provinsi Di Indonesia")
4.3. Pemetaan Distribusi Tingkat Pengangguran
ggplot(Indonesia) +
geom_sf(aes(fill = TPT22), color = "white") +
scale_fill_viridis_c(option = "plasma", name = "Tingkat Pengangguran\nTerbuka(%)") +
theme_minimal() +
labs(title = "Tingkat Pengangguran Berdasarkan\n Provinsi Di Indonesia")
5. Mencari Model Terbaik
5.1. Model Regresi OLS
modelols = lm(PPM22 ~ IPM22 + TPT22, data = Indonesia)
summary(modelols)##
## Call:
## lm(formula = PPM22 ~ IPM22 + TPT22, data = Indonesia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4388 -2.5073 0.0646 2.5480 8.2074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.1859 13.5653 5.395 6.92e-06 ***
## IPM22 -0.8578 0.2029 -4.228 0.000193 ***
## TPT22 -0.2356 0.5211 -0.452 0.654335
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.008 on 31 degrees of freedom
## Multiple R-squared: 0.4522, Adjusted R-squared: 0.4169
## F-statistic: 12.8 on 2 and 31 DF, p-value: 8.88e-05vif(modelols) #Multikolinearitas## IPM22 TPT22
## 1.287111 1.287111shapiro.test(modelols$residuals) #Normalitas Residual##
## Shapiro-Wilk normality test
##
## data: modelols$residuals
## W = 0.97344, p-value = 0.5625bptest(modelols) #Homoskedas##
## studentized Breusch-Pagan test
##
## data: modelols
## BP = 2.9191, df = 2, p-value = 0.2323AIC(modelols) #Akaike ## [1] 195.74535.2. Autokorelasi Spasial
nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=1))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=1 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 2.0555, p-value = 0.01991
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.41472476 -0.03030303 0.04687407nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=2))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=2 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.7243, p-value = 0.04233
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.23403861 -0.03030303 0.02350246nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=3))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=3 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.3565, p-value = 0.08747
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.13594311 -0.03030303 0.01502019nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=4))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=4 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.4034, p-value = 0.08026
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.11758752 -0.03030303 0.01110567nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=5))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=5 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.8396, p-value = 0.03291
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.139481565 -0.030303030 0.008517993nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=6))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=6 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.8334, p-value = 0.03338
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.121387114 -0.030303030 0.006845781nb <- knn2nb(knearneigh(st_coordinates(st_centroid(st_make_valid(Indonesia))),
k=7))
listw <- nb2listw(nb, style="W")
moran.test(resid(modelols), listw) #K=7 KNN##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: listw
##
## Moran I statistic standard deviate = 1.8317, p-value = 0.0335
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.106557639 -0.030303030 0.005582888queen= st_make_valid(Indonesia)
W=poly2nb(queen,queen=FALSE)
WL=nb2listw(W,zero.policy = TRUE)
moran.test(resid(modelols), WL) #Queen Conti##
## Moran I test under randomisation
##
## data: resid(modelols)
## weights: WL
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 1.5492, p-value = 0.06067
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.24343806 -0.03703704 0.032778305.3. Model GWR
5.3.1. GWR Fixed Gaussian
gaus=gwr.sel(PPM22 ~ TPT22+IPM22,
data=Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia)))
)
modelgaus=gwr(
PPM22 ~ TPT22+IPM22,
data = Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
bandwidth = gaus,
hatmatrix = TRUE,
se.fit = TRUE
)modelgaus$results$AICc #Akaike## [1] 192.2575prediksi <- modelgaus$SDF@data$pred
r_squared <- 1 - (sum((Indonesia$PPM22 - prediksi)^2) / sum((Indonesia$PPM22 - mean(Indonesia$PPM22))^2))
r_squared #RSquare## [1] 0.68951735.3.2. GWR Fixed Bisquare
bis=gwr.sel(PPM22 ~ TPT22+IPM22,
data=Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
gweight = gwr.bisquare,
)
modelbis=gwr(
PPM22 ~ TPT22+IPM22,
data = Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
bandwidth = bis,
hatmatrix = TRUE,
se.fit = TRUE
)modelbis$results$AICc #Akaike## [1] 193.0924prediksi <- modelbis$SDF@data$pred
r_squared <- 1 - (sum((Indonesia$PPM22 - prediksi)^2) / sum((Indonesia$PPM22 - mean(Indonesia$PPM22))^2))
r_squared #Rsquare## [1] 0.53535035.3.3. GWR Adaptive
adap=gwr.sel(PPM22 ~ TPT22+IPM22,
data=Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
adapt = TRUE)
modeladap=gwr(
PPM22 ~ TPT22+IPM22,
data = Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
adapt= adap,
hatmatrix = TRUE,
se.fit = TRUE
)modeladap$results$AICc #Akaike## [1] 253.4514prediksi <- modeladap$SDF@data$pred
r_squared <- 1 - (sum((Indonesia$PPM22 - prediksi)^2) / sum((Indonesia$PPM22 - mean(Indonesia$PPM22))^2))
r_squared #RSquare## [1] 0.83661075.3.4. GWR Fixed Tricube
tri=gwr.sel(PPM22 ~ TPT22+IPM22,
data=Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
gweight=gwr.tricube)
modeltri=gwr(
PPM22 ~ TPT22+IPM22,
data = Indonesia,
coords = st_coordinates(st_centroid(st_make_valid(Indonesia))),
bandwidth = tri,
hatmatrix = TRUE,
se.fit = TRUE
)modeltri$results$AICc #Akaike## [1] 193.3448prediksi <- modeltri$SDF@data$pred
r_squared <- 1 - (sum((Indonesia$PPM22 - prediksi)^2) / sum((Indonesia$PPM22 - mean(Indonesia$PPM22))^2))
r_squared #RSquare## [1] 0.53022766. Model Terbaik
6.1. Model GWR Fixed Gaussian
model_gwr = modelgaus
coef_gwr <- model_gwr$SDF@data[, c("(Intercept)","TPT22","IPM22")]
coef_gwr ### Koefisien Estimasi Model## (Intercept) TPT22 IPM22
## 1 39.68591 -0.34646905 -0.3935979
## 2 43.24659 -0.74544787 -0.4150339
## 3 24.16326 -0.67803371 -0.1653069
## 4 25.63020 -0.74441102 -0.1774342
## 5 25.38945 -0.72994387 -0.1752123
## 6 72.85044 -0.07164554 -0.8729795
## 7 25.82659 -0.73545994 -0.1811664
## 8 25.08100 -0.68714036 -0.1751779
## 9 26.71345 -0.74172560 -0.1923722
## 10 30.04578 -0.74123853 -0.2376397
## 11 35.64839 -0.74962982 -0.3125206
## 12 27.26369 -0.64383520 -0.2129092
## 13 39.95111 -0.66000359 -0.3813231
## 14 32.90651 -0.66646948 -0.2867009
## 15 41.49208 -0.60450491 -0.4087702
## 16 39.08067 -0.60686195 -0.3775717
## 17 23.52337 -0.60202527 -0.1634692
## 18 24.82724 -0.73009382 -0.1682307
## 19 113.68025 0.69505144 -1.4878809
## 20 103.23979 0.51465052 -1.3336889
## 21 52.39484 -0.68024511 -0.5433834
## 22 75.98104 -0.25597750 -0.8934686
## 23 131.56094 1.26637209 -1.7876902
## 24 123.10537 0.88993292 -1.6389024
## 25 26.22261 -0.62898363 -0.1948813
## 26 56.69987 -0.44271030 -0.6240666
## 27 61.90269 -0.37642778 -0.6983633
## 28 66.47480 -0.22520449 -0.7736560
## 29 73.22945 -0.14730337 -0.8685673
## 30 85.95151 0.19682220 -1.0723602
## 31 27.36320 -0.67035221 -0.2059689
## 32 24.50101 -0.70833854 -0.1659689
## 33 31.46230 -0.53625952 -0.2704359
## 34 30.93132 -0.75105738 -0.2483704model_gwr$results$AICc #Akaike## [1] 192.2575residual_gwr <- model_gwr$SDF$gwr.e
shapiro.test(residual_gwr) ##Normalitas Residual##
## Shapiro-Wilk normality test
##
## data: residual_gwr
## W = 0.98524, p-value = 0.91646.2. Plot Sebaran Untuk setiap Koefisien pada Model
6.2.1. Koefisien untuk Intercept
spplot(model_gwr$SDF, "X.Intercept.", main = "Koefisien GWR untuk Intercept", col.regions = terrain.colors(100))
6.2.2. Koefisien untuk Indeks Pembangunan Manusia
spplot(model_gwr$SDF, "IPM22", main = "Koefisien GWR untuk Indeks Pembangunan Manusia", col.regions = terrain.colors(100))
6.2.3. Koefisien untuk Tingkat Pengangguran Terbuka
spplot(model_gwr$SDF, "TPT22", main = "Koefisien GWR untuk Tingkat Pengangguran terbuka", col.regions = terrain.colors(100))
6.3. GWR Lokal dan Provinsi yang signifikan pada setiap variabel
dfree<-model_gwr$results$edf
Indonesia$IPM.t=model_gwr$SDF$IPM22/model_gwr$SDF$IPM22_se
Indonesia$IPM.t.p<-2*pt(-abs(Indonesia$IPM.t), dfree)
dfree<-model_gwr$results$edf
Indonesia$TPT.t=model_gwr$SDF$TPT22/model_gwr$SDF$TPT22_se
Indonesia$TPT.t.p<-2*pt(-abs(Indonesia$TPT.t), dfree)
# Menampilkan provinsi dengan IPM.t.p signifikan (p-value < 0.05)
prov_signif_ipm <- Indonesia[Indonesia$IPM.t.p < 0.05, ]
print(prov_signif_ipm[,"NAME_1"])## Simple feature collection with 12 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 115.821 ymin: -11.0075 xmax: 141.0194 ymax: 5.5657
## Geodetic CRS: WGS 84
## First 10 features:
## NAME_1 geometry
## 6 Gorontalo MULTIPOLYGON (((123.5491 0....
## 19 Maluku MULTIPOLYGON (((130.9183 -8...
## 20 MalukuUtara MULTIPOLYGON (((124.4479 -2...
## 21 NusaTenggaraBarat MULTIPOLYGON (((119.2006 -8...
## 22 NusaTenggaraTimur MULTIPOLYGON (((122.8815 -1...
## 23 Papua MULTIPOLYGON (((138.8911 -8...
## 24 PapuaBarat MULTIPOLYGON (((132.9052 -4...
## 26 SulawesiBarat MULTIPOLYGON (((118.996 -3....
## 27 SulawesiSelatan MULTIPOLYGON (((117.2025 -7...
## 28 SulawesiTengah MULTIPOLYGON (((123.6091 -1...# Menampilkan provinsi dengan TPT.t.p signifikan (p-value < 0.05)
prov_signif_tpt <- Indonesia[Indonesia$TPT.t.p < 0.05, ]
print(prov_signif_tpt[,"NAME_1"])## Simple feature collection with 0 features and 1 field
## Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
## Geodetic CRS: WGS 84
## [1] NAME_1 geometry
## <0 rows> (or 0-length row.names)6.4. Perbandingan Data Aktual dengan Prediksi Model Tingkat Kemiskinan
Indonesia$pred=modelgaus$SDF@data$pred
ggplot(Indonesia) +
geom_sf(aes(fill = pred), color = "black") +
scale_fill_gradientn(
colors = c("#FFF2CC", "#FFD966", "#F4B183", "#C65911"),
name = "Persentase \nPenduduk Miskin (%)",
limits = c(min(Indonesia$PPM22), max(Indonesia$PPM22))
) +
theme_minimal() +
labs(title = "Prediksi Model")+
theme(
panel.grid = element_blank(), # Remove grid lines
panel.background = element_blank(), # Remove background
axis.text = element_blank(), # Remove axis text
axis.ticks = element_blank()
)
ggplot(Indonesia) +
geom_sf(aes(fill = PPM22), color = "black") +
scale_fill_gradientn(
colors = c("#FFF2CC", "#FFD966", "#F4B183", "#C65911"),
name = "Persentase \nPenduduk Miskin (%)"
) +
theme_minimal() +
labs(title = "Data Aktual")+
theme(
panel.grid = element_blank(), # Remove grid lines
panel.background = element_blank(), # Remove background
axis.text = element_blank(), # Remove axis text
axis.ticks = element_blank()
)
7. KESIMPULAN
Metode analisis GWR dapat membantu penganalisisan datayang berupa spasial atau memiliki hubungan secarageografis. Pemodelan GWR dilakukan dengan Fixed Gaussian yang merupakan model terbaik dengan AIC sebesar 192,26 dan R2 sebesar 0,69. Pada Model OLS maupun GWR Variabel Tingkat Pegangguran Terbuka tidak signifikan terhadap model, sedangkan Variabel IPM berpengaruh signifikan pada Model OLS dan Model GWR pada beberapa Provinsi.