ANALISIS KLUSTERISASI DATA KEMISKINAN KOTA/KABUPATEN DI INDONESIA MENGGUNAKAN ALGORITMA K-MEANS PADA TAHUN 2019 - 2022
Library dan Import Data
Introduction
Dalam penelitian studi kasus ini algoritma K-Means digunakan untuk klusterisasi kota-kota/kabupaten-kabupaten di Indonesia berdasarkan variabel-variabel yang secara langsung maupun tidak langsung berhubungan dengan kemiskinan. Tujuan dari klusterisasi ini adalah untuk melihat perbandingan klusterisasi kemiskinan di kota-kota dan kabupaten-kabupaten Indonesia sebelum wabah COVID-19, awal wabah, saat wabah, dan masa akhir wabah COVID-19. Klusterisasi ini membantu dalam menentukan kota-kota/kabupaten-kabupaten yang masih termasuk ke kluster miskin dari sebelum COVID-19 sampai selesai, menentukan kota-kota/kabupaten-kabupaten yang masuk ke dalam kategori miskin karena COVID-19, dan kota-kota/kabupaten-kabupaten yang keluar dari kluster termiskin.
Tahapan Penelitian (Flowchart)
Tahapan Penelitian
Data Set
Data yang digunakan berasal dari data publikasi Badan Pusat Statistik (BPS) di tahun 2019 – 2022 yang memiliki 11 variabel yang berkaitan dengan kemiskinan.
# import data
data2019 <- read.csv2("https://raw.githubusercontent.com/juliansalomo/dataset/main/DATA2019.csv")
head(data2019)## KOTKAB AHH GK HLS IPM P0 P1 P2 PPK RLS UHH
## 1 Simeulue 63.24 404739 13.51 65.70 18.99 3.05 0.71 7210 9.08 65.22
## 2 Aceh Singkil 65.34 450217 14.30 68.91 20.78 3.83 0.98 8715 8.52 67.36
## 3 Aceh Selatan 62.31 369107 14.41 66.90 13.09 1.58 0.31 8187 8.59 64.27
## 4 Aceh Tenggara 66.00 357015 13.99 69.36 13.43 2.37 0.66 8067 9.65 68.04
## 5 Aceh Timur 66.63 417715 13.02 67.39 14.47 1.89 0.40 8600 7.86 68.67
## 6 Aceh Tengah 66.79 468577 14.26 73.14 15.50 2.08 0.46 10782 9.69 68.82
## PDRB
## 1 2211886
## 2 2395412
## 3 5479426
## 4 4906917
## 5 10280800
## 6 7472446data2020 <- read.csv2("https://raw.githubusercontent.com/juliansalomo/dataset/main/DATA2020.csv")
head(data2020)## KOTKAB AHH GK HLS IPM P0 P1 P2 PPK RLS UHH
## 1 Simeulue 63.30 444754 13.76 66.03 18.49 2.21 0.47 7085 9.34 65.26
## 2 Aceh Singkil 65.39 473983 14.31 68.94 20.20 3.55 1.03 8707 8.53 67.39
## 3 Aceh Selatan 62.40 405786 14.42 67.12 12.87 2.13 0.46 8089 8.87 64.35
## 4 Aceh Tenggara 66.11 392493 14.00 69.37 13.21 1.62 0.36 8020 9.66 68.14
## 5 Aceh Timur 66.72 440455 13.03 67.63 14.08 2.29 0.61 8489 8.15 68.72
## 6 Aceh Tengah 66.80 492227 14.27 73.24 15.08 2.89 0.95 10673 9.85 68.85
## PDRB
## 1 2274362
## 2 2422611
## 3 5530755
## 4 5058529
## 5 10605785
## 6 7387370data2021 <- read.csv2("https://raw.githubusercontent.com/juliansalomo/dataset/main/DATA2021.csv")
head(data2021)## KOTKAB AHH GK HLS IPM P0 P1 P2 PPK RLS UHH
## 1 Simeulue 63.33 458896 13.90 66.41 18.98 2.37 0.50 7148 9.48 65.28
## 2 Aceh Singkil 65.40 487249 14.32 69.22 20.36 3.67 0.91 8776 8.68 67.43
## 3 Aceh Selatan 62.46 418689 14.60 67.44 13.18 1.69 0.40 8180 8.88 64.40
## 4 Aceh Tenggara 66.17 404725 14.01 69.44 13.41 2.51 0.71 8030 9.67 68.22
## 5 Aceh Timur 66.76 460422 13.04 67.83 14.45 2.31 0.54 8577 8.21 68.74
## 6 Aceh Tengah 66.80 505933 14.28 73.37 15.26 2.88 0.80 10780 9.86 68.86
## PDRB
## 1 2440961
## 2 2702914
## 3 5960868
## 4 5401589
## 5 11654830
## 6 7999948data2022 <- read.csv2("https://raw.githubusercontent.com/juliansalomo/dataset/main/DATA2022.csv")
head(data2022)## KOTKAB AHH GK HLS IPM P0 P1 P2 PPK RLS UHH
## 1 Simeulue 63.52 493303 14.08 67.27 18.37 2.73 0.65 7371 9.73 65.48
## 2 Aceh Singkil 65.56 518951 14.34 69.62 19.18 3.18 0.74 8994 8.69 67.65
## 3 Aceh Selatan 62.68 446224 14.69 67.87 12.43 1.31 0.23 8353 8.89 64.64
## 4 Aceh Tenggara 66.39 430825 14.26 70.32 12.83 1.92 0.44 8222 9.92 68.48
## 5 Aceh Timur 66.92 491550 13.06 68.72 13.91 2.78 0.67 9127 8.32 68.94
## 6 Aceh Tengah 66.95 533810 14.61 73.95 14.50 2.41 0.59 10957 9.87 69.05
## PDRB
## 1 2650278
## 2 2984475
## 3 6426568
## 4 5813672
## 5 13008658
## 6 8827207Data Cleaning & Preparation
Pembersihan dan persiapan data adalah langkah penting dalam setiap
proyek analisis data. Bagian ini akan memastikan bahwa dataset dari 2019
hingga 2022 bebas dari nilai yang hilang missing values dan
terstruktur dengan baik, sehingga memungkinkan analisis yang akurat dan
dapat diandalkan.
#check NA
any(is.na(data2019))## [1] FALSEany(is.na(data2020))## [1] FALSEany(is.na(data2021))## [1] FALSEany(is.na(data2022))## [1] FALSE#check structure
str(data2019)## 'data.frame': 514 obs. of 12 variables:
## $ KOTKAB: chr "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ AHH : num 63.2 65.3 62.3 66 66.6 ...
## $ GK : int 404739 450217 369107 357015 417715 468577 471058 447563 462059 393198 ...
## $ HLS : num 13.5 14.3 14.4 14 13 ...
## $ IPM : num 65.7 68.9 66.9 69.4 67.4 ...
## $ P0 : num 19 20.8 13.1 13.4 14.5 ...
## $ P1 : num 3.05 3.83 1.58 2.37 1.89 2.08 3.53 2.56 2.99 1.77 ...
## $ P2 : num 0.71 0.98 0.31 0.66 0.4 0.46 0.97 0.69 0.71 0.4 ...
## $ PPK : int 7210 8715 8187 8067 8600 10782 9692 9661 9824 8889 ...
## $ RLS : num 9.08 8.52 8.59 9.65 7.86 ...
## $ UHH : num 65.2 67.4 64.3 68 68.7 ...
## $ PDRB : num 2211886 2395412 5479426 4906917 10280800 ...str(data2020)## 'data.frame': 514 obs. of 12 variables:
## $ KOTKAB: chr "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ AHH : num 63.3 65.4 62.4 66.1 66.7 ...
## $ GK : int 444754 473983 405786 392493 440455 492227 517264 477938 493145 410203 ...
## $ HLS : num 13.8 14.3 14.4 14 13 ...
## $ IPM : num 66 68.9 67.1 69.4 67.6 ...
## $ P0 : num 18.5 20.2 12.9 13.2 14.1 ...
## $ P1 : num 2.21 3.55 2.13 1.62 2.29 2.89 2.74 2.4 2.2 1.95 ...
## $ P2 : num 0.47 1.03 0.46 0.36 0.61 0.95 0.68 0.61 0.45 0.45 ...
## $ PPK : int 7085 8707 8089 8020 8489 10673 9516 9641 9816 8857 ...
## $ RLS : num 9.34 8.53 8.87 9.66 8.15 ...
## $ UHH : num 65.3 67.4 64.3 68.1 68.7 ...
## $ PDRB : num 2274362 2422611 5530755 5058529 10605785 ...str(data2021)## 'data.frame': 514 obs. of 12 variables:
## $ KOTKAB: chr "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ AHH : num 63.3 65.4 62.5 66.2 66.8 ...
## $ GK : int 458896 487249 418689 404725 460422 505933 533712 489498 504452 422685 ...
## $ HLS : num 13.9 14.3 14.6 14 13 ...
## $ IPM : num 66.4 69.2 67.4 69.4 67.8 ...
## $ P0 : num 19 20.4 13.2 13.4 14.4 ...
## $ P1 : num 2.37 3.67 1.69 2.51 2.31 2.88 3.82 2.32 3.04 2.07 ...
## $ P2 : num 0.5 0.91 0.4 0.71 0.54 0.8 1.09 0.58 0.84 0.4 ...
## $ PPK : int 7148 8776 8180 8030 8577 10780 9593 9644 9860 8867 ...
## $ RLS : num 9.48 8.68 8.88 9.67 8.21 ...
## $ UHH : num 65.3 67.4 64.4 68.2 68.7 ...
## $ PDRB : num 2440961 2702914 5960868 5401589 11654830 ...str(data2022)## 'data.frame': 514 obs. of 12 variables:
## $ KOTKAB: chr "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ AHH : num 63.5 65.6 62.7 66.4 66.9 ...
## $ GK : int 493303 518951 446224 430825 491550 533810 558638 519320 535713 451163 ...
## $ HLS : num 14.1 14.3 14.7 14.3 13.1 ...
## $ IPM : num 67.3 69.6 67.9 70.3 68.7 ...
## $ P0 : num 18.4 19.2 12.4 12.8 13.9 ...
## $ P1 : num 2.73 3.18 1.31 1.92 2.78 2.41 2.22 2.09 3.1 2.08 ...
## $ P2 : num 0.65 0.74 0.23 0.44 0.67 0.59 0.45 0.58 0.75 0.53 ...
## $ PPK : int 7371 8994 8353 8222 9127 10957 9775 9894 10211 9438 ...
## $ RLS : num 9.73 8.69 8.89 9.92 8.32 ...
## $ UHH : num 65.5 67.7 64.6 68.5 68.9 ...
## $ PDRB : num 2650278 2984475 6426568 5813672 13008658 ...Multicollinearity Test
Multikolinearitas antara variabel prediktor dapat mengganggu hasil analisis regresi. Bagian ini melakukan uji Variance Inflation Factor (VIF) pada setiap dataset untuk mengidentifikasi dan mengatasi multikolinearitas, memastikan model stabil dan dapat diinterpretasikan.
#VIF Score Test
#2019
datavif2019 <- data2019[,-1]
datavif2019$y <- c(1:nrow(data2019))
model2019<-lm(y ~., data = datavif2019)
v1<-data.frame(vif(model2019))
#2020
datavif2020<-data2020[,-1]
datavif2020$y <- c(1:nrow(data2020))
model2020<-lm(y ~., data = datavif2020)
v2<-data.frame(vif(model2020))
#2021
datavif2021<-data2021[,-1]
datavif2021$y <- c(1:nrow(data2021))
model2021<-lm(y ~., data = datavif2021)
v3<-data.frame(vif(model2021))
#2022
datavif2022<-data2022[,-1]
datavif2022$y <- c(1:nrow(data2022))
model2022<-lm(y ~., data = datavif2022)
v4<-data.frame(vif(model2022))
data.frame(cbind(v1,v2,v3,v4)) ## vif.model2019. vif.model2020. vif.model2021. vif.model2022.
## AHH 609.160911 520.429083 472.70125 392.404949
## GK 1.609971 1.596692 1.55620 1.596193
## HLS 10.949908 10.189694 10.37441 10.657364
## IPM 130.600741 125.534431 135.43524 139.041474
## P0 30.238212 30.854435 22.77435 22.860614
## P1 133.148055 142.316915 97.46023 96.164777
## P2 54.973707 58.934084 44.83959 41.478741
## PPK 20.815776 21.445255 23.07187 23.131093
## RLS 15.984060 15.391304 16.07463 16.565673
## UHH 636.277900 547.304545 498.18085 415.828081
## PDRB 1.527960 1.533444 1.55801 1.553498Dari tabel di atas, dapat dilihat variabel dengan nilai VIF tertinggi adalah Umur Harapan Hidup (UHH), sehingga variabel UHH dikeluarkan dari data.
#VIF Score Test
#2019
datavif2019<-datavif2019[,-10]
model2019<-lm(y ~., data = datavif2019)
v1<-data.frame(vif(model2019))
#2020
datavif2020<-datavif2020[,-10]
model2020<-lm(y ~., data = datavif2020)
v2<-data.frame(vif(model2020))
#2021
datavif2021<-datavif2021[,-10]
model2021<-lm(y ~., data = datavif2021)
v3<-data.frame(vif(model2021))
#2022
datavif2022<-datavif2022[,-10]
model2022<-lm(y ~., data = datavif2022)
v4<-data.frame(vif(model2022))
data.frame(cbind(v1,v2,v3,v4))## vif.model2019. vif.model2020. vif.model2021. vif.model2022.
## AHH 9.360513 9.063219 9.670079 9.971142
## GK 1.589146 1.567652 1.539841 1.581228
## HLS 10.665400 9.819720 10.042534 10.261236
## IPM 124.082402 118.438167 127.994462 130.990815
## P0 30.202091 30.722643 22.526128 22.701326
## P1 133.141738 142.116571 96.758533 95.814069
## P2 54.969879 58.841351 44.507554 41.300527
## PPK 20.317255 20.895480 22.420043 22.434112
## RLS 15.264170 14.702220 15.326874 15.767558
## PDRB 1.526607 1.533368 1.557952 1.553469Dari Tabel di atas, dapat dilihat variabel dengan nilai VIF tertinggi adalah Indeks Kedalaman Kemiskinan (P1) untuk tahun 2019 dan 2020, dan Indeks Pembangunan Manusia (IPM) untuk tahun 2021 dan 2022, sehingga variabel P1 dikeluarkan dari data 2019 dan 2020, dan IPM akan dikeluarkan dari data 2021 dan 2022. Hasil pengujian ulang nilai VIF setelah mengeluarkan variabel P1 dan IPM terdapat pada Tabel berikut.
#VIF Score Test
#2019
datavif2019<-datavif2019[,-6]
model2019<-lm(y ~., data = datavif2019)
v1<-rownames_to_column(data.frame(vif(model2019)),"var")
#2020
datavif2020<-datavif2020[,-6]
model2020<-lm(y ~., data = datavif2020)
v2<-rownames_to_column(data.frame(vif(model2020)),"var")
#2021
datavif2021<-datavif2021[,-4]
model2021<-lm(y ~., data = datavif2021)
v3<-rownames_to_column(data.frame(vif(model2021)),"var")
#2022
datavif2022<-datavif2022[,-4]
model2022<-lm(y ~., data = datavif2022)
v4<-rownames_to_column(data.frame(vif(model2022)),"var")
v1 %>% full_join(v2, by = 'var')%>% full_join(v3, by = 'var')%>% full_join(v4, by = 'var')## var vif.model2019. vif.model2020. vif.model2021. vif.model2022.
## 1 AHH 9.357730 8.963242 1.614360 1.653961
## 2 GK 1.572245 1.567039 1.495076 1.533782
## 3 HLS 10.611554 9.754322 2.894266 2.853404
## 4 IPM 123.903135 117.231890 NA NA
## 5 P0 4.419885 4.592080 20.683614 21.706475
## 6 P2 2.734130 3.043412 42.171131 40.642987
## 7 PPK 20.314715 20.472106 3.458844 3.493466
## 8 RLS 15.261215 14.589985 3.887971 3.907663
## 9 PDRB 1.526217 1.533208 1.445400 1.450150
## 10 P1 NA NA 92.953319 94.891084Dari Tabel di atas, dapat dilihat variabel dengan nilai VIF tertinggi adalah Indeks Pembangunan Manusia (IPM) untuk tahun 2019 dan 2020, dan Indeks Kedalaman Kemiskinan (P1) untuk tahun 2021 dan 2022, sehingga variabel IPM dikeluarkan dari data 2019 dan 2020, dan P1 dikeluarkan dari data 2021 dan 2022. Hasil pengujian ulang nilai VIF setelah mengeluarkan variabel P1 dan IPM terdapat pada Tabel di bawah.
#VIF Score Test
#2019
datavif2019<-datavif2019[,-4]
model2019<-lm(y ~., data = datavif2019)
v1<-rownames_to_column(data.frame(vif(model2019)),"var")
#2020
datavif2020<-datavif2020[,-4]
model2020<-lm(y ~., data = datavif2020)
v2<-rownames_to_column(data.frame(vif(model2020)),"var")
#2021
datavif2021<-datavif2021[,-5]
model2021<-lm(y ~., data = datavif2021)
v3<-rownames_to_column(data.frame(vif(model2021)),"var")
#2022
datavif2022<-datavif2022[,-5]
model2022<-lm(y ~., data = datavif2022)
v4<-rownames_to_column(data.frame(vif(model2022)),"var")
v1 %>% full_join(v2, by = 'var')%>% full_join(v3, by = 'var')%>% full_join(v4, by = 'var')## var vif.model2019. vif.model2020. vif.model2021. vif.model2022.
## 1 AHH 1.646312 1.611087 1.614227 1.652693
## 2 GK 1.542270 1.522620 1.480975 1.532910
## 3 HLS 2.898568 2.867352 2.894253 2.847064
## 4 P0 4.053101 4.277329 3.635735 3.968297
## 5 P2 2.713438 3.025343 2.299201 2.687676
## 6 PPK 3.529611 3.338804 3.441101 3.475850
## 7 RLS 3.837108 3.777963 3.836208 3.871524
## 8 PDRB 1.431771 1.432858 1.444330 1.450074Dari Tabel di atas, semua VIF untuk tiap tahun sudah di bawah 10, sehingga gejala multikolinieritas untuk data penelitian telah diatasi, dan sudah tidak ada gejala multikolinieritas lagi.
Data Normalization
Sebelum melanjutkan proses klusterisasi, peneliti melakukan transformasi data menggunakan PCA untuk mengekstrak komponen-komponen penting atau principal component untuk mengurangi variabel yang tidak perlu.
Sebelum menerapkan PCA, peneliti perlu melakukan normalisasi data terlebih dahulu.
# Data after Removing Variables Detected with Multicollinearity
datascaled2019 <- datavif2019[,-9] %>% mutate_at(c(1:ncol(datavif2019[,-9])), list(~c(scale(.))))
datascaled2020 <- datavif2020[,-9] %>% mutate_at(c(1:ncol(datavif2020[,-9])), list(~c(scale(.))))
datascaled2021 <- datavif2021[,-9] %>% mutate_at(c(1:ncol(datavif2021[,-9])), list(~c(scale(.))))
datascaled2022 <- datavif2022[,-9] %>% mutate_at(c(1:ncol(datavif2022[,-9])), list(~c(scale(.))))Selain data tanpa gejala multikolinieritas, peneliti juga menggunakan data yang masih memiliki gejala multikolinieritas untuk membandingkan hasil klusterisasi yang lebih baik.
# Data with Multicollinearity Symptoms
datascaled2019X <- data2019[,-1] %>% mutate_at(c(1:ncol(data2019[,-1])), list(~c(scale(.))))
datascaled2020X <- data2020[,-1] %>% mutate_at(c(1:ncol(data2020[,-1])), list(~c(scale(.))))
datascaled2021X <- data2021[,-1] %>% mutate_at(c(1:ncol(data2021[,-1])), list(~c(scale(.))))
datascaled2022X <- data2022[,-1] %>% mutate_at(c(1:ncol(data2022[,-1])), list(~c(scale(.))))PCA Principal Component Analysis
Setelah dinormalisasi peneliti melanjutkan pemilihan jumlah komponen yang akan diambil atau diekstrak berdasarkan Tabel di bawah
# Data without Multicollinearity Symptoms
datapca2019 <- prcomp(datascaled2019)
datapca2020 <- prcomp(datascaled2020)
datapca2021 <- prcomp(datascaled2021)
datapca2022 <- prcomp(datascaled2022)
# Data with Multicollinearity Symptoms
datapca2019X <- prcomp(datascaled2019X)
datapca2020X <- prcomp(datascaled2020X)
datapca2021X <- prcomp(datascaled2021X)
datapca2022X <- prcomp(datascaled2022X)a <- cbind(Tahun=2019,summary(datapca2019)$importance%>%data.frame())
b <- cbind(Tahun=2020,summary(datapca2020)$importance%>%data.frame())
c <- cbind(Tahun=2021,summary(datapca2021)$importance%>%data.frame())
d <- cbind(Tahun=2022,summary(datapca2022)$importance%>%data.frame())
rbind(a,b,c,d)## Tahun PC1 PC2 PC3 PC4 PC5
## Standard deviation 2019 1.976898 1.179840 0.9783243 0.7763630 0.7314797
## Proportion of Variance 2019 0.488520 0.174000 0.1196400 0.0753400 0.0668800
## Cumulative Proportion 2019 0.488520 0.662520 0.7821600 0.8575000 0.9243800
## Standard deviation1 2020 1.967209 1.185989 0.9859026 0.7738675 0.7485962
## Proportion of Variance1 2020 0.483740 0.175820 0.1215000 0.0748600 0.0700500
## Cumulative Proportion1 2020 0.483740 0.659560 0.7810600 0.8559200 0.9259700
## Standard deviation2 2021 1.958563 1.160546 0.9869264 0.7932754 0.7518250
## Proportion of Variance2 2021 0.479500 0.168360 0.1217500 0.0786600 0.0706600
## Cumulative Proportion2 2021 0.479500 0.647850 0.7696100 0.8482700 0.9189200
## Standard deviation3 2022 1.965557 1.186083 0.9912827 0.7738107 0.7458124
## Proportion of Variance3 2022 0.482930 0.175850 0.1228300 0.0748500 0.0695300
## Cumulative Proportion3 2022 0.482930 0.658780 0.7816100 0.8564500 0.9259800
## PC6 PC7 PC8
## Standard deviation 0.5323218 0.4213806 0.3794789
## Proportion of Variance 0.0354200 0.0222000 0.0180000
## Cumulative Proportion 0.9598000 0.9820000 1.0000000
## Standard deviation1 0.5176301 0.4341814 0.3685019
## Proportion of Variance1 0.0334900 0.0235600 0.0169700
## Cumulative Proportion1 0.9594600 0.9830300 1.0000000
## Standard deviation2 0.5545195 0.4433964 0.3801603
## Proportion of Variance2 0.0384400 0.0245800 0.0180700
## Cumulative Proportion2 0.9573600 0.9819300 1.0000000
## Standard deviation3 0.5068357 0.4353078 0.3817790
## Proportion of Variance3 0.0321100 0.0236900 0.0182200
## Cumulative Proportion3 0.9580900 0.9817800 1.0000000Melihat dari proportion of variance pada Tabel diatas, Principal Component Analysis, peneliti memilih 2 komponen pertama, yaitu PC1 dan PC2, karena peneliti melihat bahwa 2 komponen tersebut merupakan komponen paling penting dengan proportion of variance terbesar.
Selain itu peneliti juga melihat dari eigenvalue untuk menentukan jumlah komponen yang akan dipilih mengikuti Kaiser’s rule. Adapun eigenvalue untuk data PCA di atas terdapat pada Tabel berikut.
data.frame(ev2019 = get_eigenvalue(datapca2019)$eigenvalue,
ev2020 = get_eigenvalue(datapca2020)$eigenvalue,
ev2021 = get_eigenvalue(datapca2021)$eigenvalue,
ev2022 = get_eigenvalue(datapca2022)$eigenvalue) %>% t() %>% data.frame() %>% round(2)## X1 X2 X3 X4 X5 X6 X7 X8
## ev2019 3.91 1.39 0.96 0.60 0.54 0.28 0.18 0.14
## ev2020 3.87 1.41 0.97 0.60 0.56 0.27 0.19 0.14
## ev2021 3.84 1.35 0.97 0.63 0.57 0.31 0.20 0.14
## ev2022 3.86 1.41 0.98 0.60 0.56 0.26 0.19 0.15Dari Tabel di atas, eigenvalue yang lebih besar dari 1 untuk setiap tahunnya adalah komponen pertama (PC1) dan komponen kedua (PC2), sehingga peneliti memilih 2 komponen pertama.
Selain data tanpa gejala multikolinieritas, peneliti juga menggunakan data yang masih memiliki gejala multikolinieritas untuk membandingkan hasil klusterisasi yang lebih baik. Adapun eigenvalue untuk data PCA dari data yang masih ada gejala multikolinieritas terdapat pada Tabel berikut.
a <- cbind(Tahun=2019,summary(datapca2019X)$importance%>%data.frame())
b <- cbind(Tahun=2020,summary(datapca2020X)$importance%>%data.frame())
c <- cbind(Tahun=2021,summary(datapca2021X)$importance%>%data.frame())
d <- cbind(Tahun=2022,summary(datapca2022X)$importance%>%data.frame())
rbind(a,b,c,d)## Tahun PC1 PC2 PC3 PC4 PC5
## Standard deviation 2019 2.451161 1.298372 1.077297 0.9398396 0.7897986
## Proportion of Variance 2019 0.546200 0.153250 0.105510 0.0803000 0.0567100
## Cumulative Proportion 2019 0.546200 0.699450 0.804960 0.8852600 0.9419600
## Standard deviation1 2020 2.439933 1.307624 1.080959 0.9540030 0.7945264
## Proportion of Variance1 2020 0.541210 0.155440 0.106220 0.0827400 0.0573900
## Cumulative Proportion1 2020 0.541210 0.696650 0.802880 0.8856100 0.9430000
## Standard deviation2 2021 2.430418 1.278533 1.083126 0.9668936 0.8150851
## Proportion of Variance2 2021 0.536990 0.148600 0.106650 0.0849900 0.0604000
## Cumulative Proportion2 2021 0.536990 0.685600 0.792250 0.8772400 0.9376400
## Standard deviation3 2022 2.441285 1.303753 1.081919 0.9526824 0.7957976
## Proportion of Variance3 2022 0.541810 0.154520 0.106410 0.0825100 0.0575700
## Cumulative Proportion3 2022 0.541810 0.696330 0.802740 0.8852500 0.9428300
## PC6 PC7 PC8 PC9 PC10
## Standard deviation 0.5497785 0.4218793 0.3829174 0.07773412 0.06839940
## Proportion of Variance 0.0274800 0.0161800 0.0133300 0.00055000 0.00043000
## Cumulative Proportion 0.9694400 0.9856200 0.9989500 0.99950000 0.99993000
## Standard deviation1 0.5374245 0.4346648 0.3708697 0.07981078 0.06612419
## Proportion of Variance1 0.0262600 0.0171800 0.0125000 0.00058000 0.00040000
## Cumulative Proportion1 0.9692600 0.9864300 0.9989400 0.99952000 0.99991000
## Standard deviation2 0.5709770 0.4456901 0.3842727 0.08784789 0.07032233
## Proportion of Variance2 0.0296400 0.0180600 0.0134200 0.00070000 0.00045000
## Cumulative Proportion2 0.9672700 0.9853300 0.9987600 0.99946000 0.99991000
## Standard deviation3 0.5262161 0.4361773 0.3849776 0.08370828 0.07283565
## Proportion of Variance3 0.0251700 0.0173000 0.0134700 0.00064000 0.00048000
## Cumulative Proportion3 0.9680000 0.9853000 0.9987700 0.99941000 0.99989000
## PC11
## Standard deviation 0.02833024
## Proportion of Variance 0.00007000
## Cumulative Proportion 1.00000000
## Standard deviation1 0.03057783
## Proportion of Variance1 0.00009000
## Cumulative Proportion1 1.00000000
## Standard deviation2 0.03205606
## Proportion of Variance2 0.00009000
## Cumulative Proportion2 1.00000000
## Standard deviation3 0.03513288
## Proportion of Variance3 0.00011000
## Cumulative Proportion3 1.00000000data.frame(ev2019 = get_eigenvalue(datapca2019X)$eigenvalue,
ev2020 = get_eigenvalue(datapca2020X)$eigenvalue,
ev2021 = get_eigenvalue(datapca2021X)$eigenvalue,
ev2022 = get_eigenvalue(datapca2022X)$eigenvalue) %>% t() %>% data.frame() %>% round(2)## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11
## ev2019 6.01 1.69 1.16 0.88 0.62 0.30 0.18 0.15 0.01 0.00 0
## ev2020 5.95 1.71 1.17 0.91 0.63 0.29 0.19 0.14 0.01 0.00 0
## ev2021 5.91 1.63 1.17 0.93 0.66 0.33 0.20 0.15 0.01 0.00 0
## ev2022 5.96 1.70 1.17 0.91 0.63 0.28 0.19 0.15 0.01 0.01 0Dari Tabel di atas, eigenvalue yang lebih besar dari 1 untuk setiap tahunnya adalah komponen pertama (PC1) sampai komponen ketiga (PC3), sehingga peneliti memilih 3 komponen pertama.
#Optimizing PCA Selection and Data Finalization for Data without Multicollinearity Symptoms
datafinal2019 <- data.frame(datapca2019$x[,1:2])
datafinal2020 <- data.frame(datapca2020$x[,1:2])
datafinal2021 <- data.frame(datapca2021$x[,1:2])
datafinal2022 <- data.frame(datapca2022$x[,1:2])
rownames(datafinal2019) <- data2019$KOTKAB
rownames(datafinal2020) <- data2020$KOTKAB
rownames(datafinal2021) <- data2021$KOTKAB
rownames(datafinal2022) <- data2022$KOTKAB
#Optimizing PCA Selection and Data Finalization for Data with Multicollinearity Symptoms
datafinal2019X <- data.frame(datapca2019X$x[,1:3])
datafinal2020X <- data.frame(datapca2020X$x[,1:3])
datafinal2021X <- data.frame(datapca2021X$x[,1:3])
datafinal2022X <- data.frame(datapca2022X$x[,1:3])
rownames(datafinal2019X) <- data2019$KOTKAB
rownames(datafinal2020X) <- data2020$KOTKAB
rownames(datafinal2021X) <- data2021$KOTKAB
rownames(datafinal2022X) <- data2022$KOTKABData Clustering
Data Optimum Cluster Selection
Untuk penentuan jumlah kluster peneliti menggunakan bantuan 3 metode yang dapat membantu penentuan jumlah kluster, yaitu Metode Elbow, Metode Silhouette, dan Metode Gap Statistic
Elbow Method
Hasil jumlah kluster optimal data tanpa gejala multikolinieritas, berdasarkan 8 variabel menggunakan metode elbow dapat dilihat pada Gambar berikut
#Elbow Method for Cluster Number Selection in Data without Multicollinearity Symptoms
e1 <- fviz_nbclust(datafinal2019, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2019")
e2 <- fviz_nbclust(datafinal2020, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2020")
e3 <- fviz_nbclust(datafinal2021, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2021")
e4 <- fviz_nbclust(datafinal2022, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2022")
grid.arrange(e1, e2, e3, e4, nrow = 2)
Elbow Method for Cluster Number Selection in Data without Multicollinearity Symptoms
Dari Gambar di atas, terlihat adanya lekukan tepat setelah kluster 3 yang membentuk bentuk elbow atau siku, sehingga berdasarkan metode elbow jumlah kluster optimal untuk data tanpa gejala multikolinieritas adalah 3 kluster.
Selanjutnya untuk data dengan gejala multikolinieritas, berdasarkan 11 variabel yang ada, hasil jumlah kluster optimal menggunakan metode elbow dapat dilihat pada Gambar berikut.
#Elbow Method for Cluster Number Selection in Data with Multicollinearity Symptoms
e1x <- fviz_nbclust(datafinal2019X, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2019")
e2x <- fviz_nbclust(datafinal2020X, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2020")
e3x <- fviz_nbclust(datafinal2021X, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2021")
e4x <- fviz_nbclust(datafinal2022X, kmeans, method = "wss") +
labs(subtitle = "Elbow Method 2022")
grid.arrange(e1x, e2x, e3x, e4x, nrow = 2)
Elbow Method for Cluster Number Selection in Data with Multicollinearity Symptoms
Dari Gambar di atas, terlihat adanya lekukan tepat setelah kluster 3 yang membentuk bentuk elbow atau siku, sehingga berdasarkan metode elbow jumlah kluster optimal untuk data dengan gejala multikolinieritas adalah 3 kluster.
Silhouette Method
Hasil jumlah kluster optimal untuk data tanpa gejala multikolinieritas menggunakan metode Silhouette dapat dilihat pada Gambar berikut.
#Silhouette Method for Cluster Number Selection in Data without Multicollinearity Symptoms
s1<- fviz_nbclust(datafinal2019, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2019")
s2<-fviz_nbclust(datafinal2020, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2020")
s3<-fviz_nbclust(datafinal2021, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2021")
s4 <- fviz_nbclust(datafinal2022, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2022")
grid.arrange(s1, s2, s3, s4, nrow = 2)
Silhouette Method for Cluster Number Selection in Data without Multicollinearity Symptoms
Dari Gambar di atas, terlihat bahwa jumlah kluster optimal untuk tahun 2019 sampai 2021 adalah 3 kluster, sedangkan untuk tahun 2022 adalah 2 kluster. Tapi karena di sini peneliti ingin membandingkan hasil kluster dari keempat tahun, maka peneliti memilih jumlah kluster yang sama. Dan berdasarkan Gambar di atas, jumlah kluster 2 dan 3 di tahun 2022 memiliki nilai indeks silhouette yang tidak berbeda jauh, sehingga berdasarkan metode silhouette jumlah kluster optimal untuk data tanpa gejala multikolinieritas adalah 3 kluster.
Selanjutnya untuk data dengan gejala multikolinieritas, berdasarkan 11 variabel yang ada, hasil jumlah kluster optimal menggunakan metode silhouette dapat dilihat pada Gambar berikut.
#Silhouette Method for Cluster Number Selection in Data with Multicollinearity Symptoms
s1X<- fviz_nbclust(datafinal2019X, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2019")
s2X<-fviz_nbclust(datafinal2020X, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2020")
s3X<-fviz_nbclust(datafinal2021X, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2021")
s4X <- fviz_nbclust(datafinal2022X, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette Method 2022")
grid.arrange(s1X, s2X, s3X, s4X, nrow = 2)
Silhouette Method for Cluster Number Selection in Data with Multicollinearity Symptoms
Dari Gambar di atas, terlihat bahwa jumlah kluster optimal untuk tahun 2019 sampai 2021 adalah 3 kluster, sedangkan untuk tahun 2022 adalah 2 kluster. Tapi karena di sini peneliti ingin membandingkan hasil kluster dari keempat tahun, maka peneliti memilih jumlah kluster yang sama. Dan berdasarkan Gambar di atas jumlah kluster 2 dan 3 di tahun 2022 memiliki nilai indeks silhouette yang tidak berbeda jauh, sehingga berdasarkan metode silhouette jumlah kluster optimal untuk data dengan gejala multikolinieritas adalah 3 kluster.
Gap Statistic Method
Hasil jumlah kluster optimal untuk data tanpa gejala multikolinieritas menggunakan metode Gap Statistic dapat dilihat pada Gambar berikut.
#Gap Statistic Method for Cluster Number Selection in Data without Multicollinearity Symptoms
gs1 <- fviz_nbclust(datafinal2019, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2019")
gs2 <- fviz_nbclust(datafinal2020, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2020")
gs3 <- fviz_nbclust(datafinal2021, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2021")
gs4 <- fviz_nbclust(datafinal2022, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2022")
grid.arrange(gs1, gs2, gs3, gs4, nrow = 2)
Gap Statistic Method for Cluster Number Selection in Data without Multicollinearity Symptoms
Dari Gambar di atas, terlihat bahwa jumlah kluster optimal untuk tahun 2019 sampai tahun 2022 adalah 1 kluster, sehingga berdasarkan metode gap statistic jumlah kluster optimal untuk data tanpa gejala multikolinieritas adalah 1 kluster.
Selanjutnya untuk data dengan gejala multikolinieritas, berdasarkan 11 variabel yang ada, hasil jumlah kluster optimal menggunakan metode gap statistic dapat dilihat pada Gambar berikut.
#Gap Statistic Method for Cluster Number Selection in Data with Multicollinearity Symptoms
gs1X <- fviz_nbclust(datafinal2019X, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2019")
gs2X <- fviz_nbclust(datafinal2020X, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2020")
gs3X <- fviz_nbclust(datafinal2021X, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2021")
gs4X <- fviz_nbclust(datafinal2022X, kmeans, method = "gap_stat")+
labs(subtitle = "Gap Statistic Method 2022")## Warning: did not converge in 10 iterationsgrid.arrange(gs1X, gs2X, gs3X, gs4X, nrow = 2)
Gap Statistic Method for Cluster Number Selection in Data with Multicollinearity Symptoms
Dari Gambar di atas, terlihat bahwa jumlah kluster optimal untuk tahun 2019 sampai tahun 2022 adalah 1 kluster, sehingga berdasarkan metode gap statistic jumlah kluster optimal untuk data dengan gejala multikolinieritas adalah 1 kluster.
Karena Elbow dan Silhoutte menunjukkan jumlah kluster optimum adalah 3, dan hanya Gap Statistic yang memiliki hasil berbeda, Maka kita memutuskan jumlah kluster optimum adalah 3
K-Means Clustering
Untuk mengetahui seberapa baik hasil klusterisasi k-means terhadap data yang digunakan, maka perlu memeriksa plot silhouette dari hasil klusterisasi. Gambar berikut adalah plot silhouette untuk data tanpa gejala multikolinieritas.
# K-Means Clustering in data without Multicollinearity Symptoms
km2019 <- eclust(datafinal2019, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2020 <- eclust(datafinal2020, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2021 <- eclust(datafinal2021, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2022 <- eclust(datafinal2022, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2019$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2019$cluster)]%>%as.integer()
km2020$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2020$cluster)]%>%as.integer()
km2021$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2021$cluster)]%>%as.integer()
km2022$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2022$cluster)]%>%as.integer()
grid.arrange(fviz_silhouette(km2019, print.summary = F)+ labs(subtitle = "2019: 3 Clusters without Multicollinearity Symptoms"),
fviz_silhouette(km2020, print.summary = F)+ labs(subtitle = "2020: 3 Clusters without Multicollinearity Symptoms"),
fviz_silhouette(km2021, print.summary = F)+ labs(subtitle = "2021: 3 Clusters without Multicollinearity Symptoms"),
fviz_silhouette(km2022, print.summary = F)+ labs(subtitle = "2022: 3 Clusters without Multicollinearity Symptoms"),nrow = 2)
K-Means Clustering in data without Multicollinearity Symptoms
Dari Gambar di atas, dapat dilihat kluster di tahun 2019 sampai 2022 memiliki nilai di atas nilai rata-rata, lalu hanya ada sedikit dari hasil kluster yang memiliki nilai negatif. Selain itu nilai rata-rata dari setiap tahun berada di atas 0,5, yang mana jika dibulatkan lebih mendekati nilai 1 dibanding 0, sehingga dapat dikatakan data terbagi dengan baik.
Selanjutnya Gambar berikut adalah plot silhouette untuk data dengan gejala multikolinieritas.
# K-Means Clustering in data with Multicollinearity Symptoms
km2019X <- eclust(datafinal2019X, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2020X <- eclust(datafinal2020X, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2021X <- eclust(datafinal2021X, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2022X <- eclust(datafinal2022X, k=3, FUNcluster="kmeans", hc_metric="euclidean", graph=F)
km2019X$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2019X$cluster)]%>%as.integer()
km2020X$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2020X$cluster)]%>%as.integer()
km2021X$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2021X$cluster)]%>%as.integer()
km2022X$cluster <- c("1"="2", "2" = "3", "3" = "1")[as.factor(km2022X$cluster)]%>%as.integer()
grid.arrange(fviz_silhouette(km2019X, print.summary = F)+ labs(subtitle = "2019: 3 Clusters with Multicollinearity Symptoms"),
fviz_silhouette(km2020X, print.summary = F)+ labs(subtitle = "2020: 3 Clusters with Multicollinearity Symptoms"),
fviz_silhouette(km2021X, print.summary = F)+ labs(subtitle = "2021: 3 Clusters with Multicollinearity Symptoms"),
fviz_silhouette(km2022X, print.summary = F)+ labs(subtitle = "2022: 3 Clusters with Multicollinearity Symptoms"),nrow = 2)
K-Means Clustering in data with Multicollinearity Symptoms
Dari Gambar di atas, dapat dilihat kluster di tahun 2019 sampai 2022 memiliki nilai di atas nilai rata-rata, lalu hanya ada sedikit dari hasil kluster yang memiliki nilai negatif. Namun nilai rata-rata dari setiap tahun berada di bawah 0,5, yang mana jika dibulatkan lebih mendekati nilai 0 dibanding 1, sehingga dapat dikatakan data terbagi kurang baik.
Menggunakan algoritma K-Means dengan jumlah kluster 3, Gambar berikut adalah plot hasil kluster dari data tanpa gejala multikolinieritas.
# Cluster Plot in data with Multicollinearity Symptoms
clusterkm2019 <- fviz_cluster(km2019, data = datafinal2019, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2019")
clusterkm2020 <- fviz_cluster(km2020, data = datafinal2020, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2020")
clusterkm2021 <- fviz_cluster(km2021, data = datafinal2021, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2021")
clusterkm2022 <- fviz_cluster(km2022, data = datafinal2022, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2022")
plot<-grid.arrange(clusterkm2019, clusterkm2020, clusterkm2021, clusterkm2022, nrow=2)
Cluster Plot in data without Multicollinearity Symptoms
Dari plot kluster di Gambar di atas, dapat dilihat bahwa tiap kluster terbagi dengan jelas. Selanjutnya plot hasil kluster untuk data dengan gejala multikolinieritas dapat dilihat pada Gambar berikut.
# Cluster Plot in data with Multicollinearity Symptoms
clusterkm2019 <- fviz_cluster(km2019X, data = datafinal2019, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2019")
clusterkm2020 <- fviz_cluster(km2020X, data = datafinal2020, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2020")
clusterkm2021 <- fviz_cluster(km2021X, data = datafinal2021, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2021")
clusterkm2022 <- fviz_cluster(km2022X, data = datafinal2022, labelsize = 0) +
labs(subtitle = "K-Means With PCA & K-Value = 3 in 2022")
grid.arrange(clusterkm2019, clusterkm2020, clusterkm2021, clusterkm2022, nrow=2)
Cluster Plot in data with Multicollinearity Symptoms
Dari plot kluster di Gambar 4.10, dapat dilihat bahwa tiap kluster terbagi dengan buruk. Di tahun 2019 dan 2021 hasil kluster ketiga kluster saling bertumpukan, di tahun 2020 dan 2022 ada irisan antara kluster 1 dengan 2 dan 2 dengan 2. Maka dengan ini data dengan gejala multikolinieritas tidak akan dianalisis lebih lanjut, dan penelitian akan fokus dengan data tanpa gejala multikolinieritas.
Data Profile
Tabel berikut adalah ringkasan data hasil klusterisasi.
#Assign into new object
dataprofile2019 <- data2019
dataprofile2020 <- data2020
dataprofile2021 <- data2021
dataprofile2022 <- data2022
#Assign cluster result into the new object
dataprofile2019$cluster <- km2019$cluster
dataprofile2020$cluster <- km2020$cluster
dataprofile2021$cluster <- km2021$cluster
dataprofile2022$cluster <- km2022$cluster
#Print cluster summary
rbind(data.frame(Tahun=2019,dataprofile2019[,-1] %>%
group_by(cluster) %>%
summarise_all(mean)),
data.frame(Tahun=2020,dataprofile2020[,-1] %>%
group_by(cluster) %>%
summarise_all(mean)),
data.frame(Tahun=2021,dataprofile2021[,-1] %>%
group_by(cluster) %>%
summarise_all(mean)),
data.frame(Tahun=2022,dataprofile2022[,-1] %>%
group_by(cluster) %>%
summarise_all(mean))) %>% round(2)## Tahun cluster AHH GK HLS IPM P0 P1 P2 PPK RLS
## 1 2019 1 70.57 512287.9 14.06 77.70 6.55 0.98 0.24 13871.60 10.20
## 2 2019 2 66.95 376294.9 12.75 68.30 11.63 1.78 0.43 9769.85 7.86
## 3 2019 3 62.21 475376.5 10.55 56.03 32.55 8.19 2.98 6413.33 5.59
## 4 2020 1 70.83 550931.2 14.20 78.13 6.43 0.96 0.24 13820.97 10.40
## 5 2020 2 67.18 405488.0 12.83 68.58 11.65 1.77 0.42 9665.83 8.03
## 6 2020 3 62.35 507387.5 10.57 55.62 31.99 7.89 2.77 6078.18 5.63
## 7 2021 1 70.99 577015.6 14.29 78.60 6.80 1.00 0.25 13975.57 10.57
## 8 2021 2 67.38 423560.3 12.90 69.00 11.69 1.83 0.45 9790.05 8.12
## 9 2021 3 62.51 512906.7 10.97 56.98 31.34 6.99 2.41 6216.25 6.06
## 10 2022 1 71.14 610589.2 14.31 79.06 6.43 0.99 0.25 14296.97 10.65
## 11 2022 2 67.58 449343.0 12.97 69.61 11.08 1.71 0.42 10085.78 8.22
## 12 2022 3 62.88 553122.9 11.07 57.86 30.69 6.91 2.35 6427.27 6.22
## UHH PDRB
## 1 72.51 82901971
## 2 68.94 17972485
## 3 64.15 4127070
## 4 72.73 86586677
## 5 69.16 18122001
## 6 64.26 3450694
## 7 72.88 92567657
## 8 69.36 20432609
## 9 64.42 3695297
## 10 73.07 101592981
## 11 69.60 23112767
## 12 64.86 3978529Dari Tabel di atas, dapat dilihat bahwa kluster pertama merupakan kluster dengan tingkat kemiskinan terendah karena memiliki nilai tertinggi dalam pendapatan domestik regional bruto, pendapatan per kapita, indeks pembangunan manusia, garis kemiskinan, angka harapan hidup, umur harapan hidup, harapan lama sekolah dan rata-rata lama sekolah dan memiliki nilai terendah dalam persentase penduduk miskin, kedalaman kemiskinan dan keparahan kemiskinan.
Sebaliknya, kluster ketiga merupakan kluster dengan tingkat kemiskinan tertinggi karena memiliki nilai kebalikan dari kluster pertama, sehingga kluster kedua merupakan kluster dengan tingkat kemiskinan sedang.
Untuk ringkasan jumlah kota/kabupaten tiap kluster pada tahun 2019 sampai 2022 terdapat pada Tabel berikut.
c2019 <- dataprofile2019 %>% count(cluster, name = 'c2019')
c2020 <- dataprofile2020 %>% count(cluster, name = 'c2020')
c2021 <- dataprofile2021 %>% count(cluster, name = 'c2021')
c2022 <- dataprofile2022 %>% count(cluster, name = 'c2022')
df<-c2019 %>% full_join(c2020, by = 'cluster') %>% full_join(c2021, by = 'cluster') %>% full_join(c2022, by = 'cluster')
df%>%
pivot_longer(cols=c(-cluster),names_to="Year")%>%
pivot_wider(names_from=c(cluster))## # A tibble: 4 × 4
## Year `1` `2` `3`
## <chr> <int> <int> <int>
## 1 c2019 114 364 36
## 2 c2020 103 377 34
## 3 c2021 100 374 40
## 4 c2022 103 371 40Dari Tabel 4.17, dapat dilihat bahwa dari tahun 2019 ke tahun 2020 adanya penurunan kota/kabupaten di kluster 1 dan kluster 3, sebaliknya ada penambahan di kluster 2 yang menandakan adanya kota/kabupaten yang mengalami perkembangan dan ada yang mengalami penurunan. Selanjutnya dari tahun 2020 ke tahun 2021, terjadi penurunan jumlah kota/kabupaten di kluster 1 dan 2, dan ada penambahan jumlah kluster di kluster 3 yang menandakan dari tahun 2020 ke 2021 diduga beberapa kota/kabupaten terkena dampak dari wabah COVID-19 yang menyebabkan mereka mengelami penurunan dan berpindah kluster dari kluster 1 ke kluster 2, dari kluster 2 ke kluster 3. Dan yang terakhir dari tahun 2021 ke tahun 2022, terjadi penambahan di kluster 1 dan pengurangan di kluster 2 yang menandakan adanya kota/kabupaten yang mengalami perkembangan. Jika dilihat dari awal periode ke akhir periode, maka dapat dilihat bahwa ada kota/kabupaten yang mengalami penurunan, karena kluster 1, kluster dengan tingkat kemiskinan terendah yang awalnya berisi 114 kota/kabupaten, menurun menjadi 103 kota/kabupaten. Sebaliknya kluster 2 dan 3 mengalami kenaikan dari 364 dan 36 kota/kabupaten menjadi 371 dan 40 kota/kabupaten.
Selanjutnya peneliti meneliti kota/kabupaten yang mengalami penurunan dan yang mengalami perkembangan di masa COVID-19. Adapun kota/kabupaten yang mengalami penurunan di akhir masa COVID-19 (2022) berdasarkan masa sebelum COVID-19 (2019) terdapat pada Tabel berikut.
cluster_2019_to_2022 <- data.frame(KOTKAB = data2019$KOTKAB,
c2019 = dataprofile2019$cluster,
c2020 = dataprofile2020$cluster,
c2021 = dataprofile2021$cluster,
c2022 = dataprofile2022$cluster)
#Mengalami penurunan di akhir masa CoVid-19 berdasarkan masa sebelum CoVid-19
cluster_2019_to_2022%>% filter(c2019<c2022) %>% arrange(c2022)## KOTKAB c2019 c2020 c2021 c2022
## 1 Kampar 1 2 2 2
## 2 Kota Lubuklinggau 1 2 2 2
## 3 Bangka 1 1 2 2
## 4 Kep. Seribu 1 2 2 2
## 5 Bogor 1 1 2 2
## 6 Kota Cirebon 1 2 2 2
## 7 Klaten 1 2 2 2
## 8 Karanganyar 1 2 2 2
## 9 Kota Pekalongan 1 2 2 2
## 10 Kota Probolinggo 1 2 2 2
## 11 Jayapura 1 2 2 2
## 12 Sumba Barat 2 2 3 3
## 13 Lembata 2 2 3 3
## 14 Rote Ndao 2 2 3 3
## 15 Manggarai Timur 2 2 3 3
## 16 Maluku Tenggara Barat 2 3 3 3Dari Tabel di atas, dapat dilihat ada 16 kota/kabupaten yang mengalami penurunan di akhir masa COVID-19 berdasarkan masa sebelum COVID-19. 11 kota/kabupaten mengalami penurunan dari kluster kemiskinan terendah menjadi kluster dengan tingkat kemiskinan sedang, sedangkan 5 kota/kabupaten mengalami penurunan dari kluster kemiskinan tingkat sedang menjadi kluster dengan tingkat kemiskinan tertinggi.
Selanjutnya, kota/kabupaten yang mengalami perkembangan di akhir masa COVID-19 (2022) berdasarkan masa sebelum COVID-19 (2019) terdapat pada Tabel berikut.
#Mengalami perkembangan di akhir masa CoVid-19 berdasarkan masa sebelum CoVid-19
cluster_2019_to_2022%>% filter(c2019>c2022) %>% arrange(c2022)## KOTKAB c2019 c2020 c2021 c2022
## 1 Kepulauan Meranti 3 2 2 2Dari Tabel di atas dapat dilihat bahwa hanya ada 1 kabupaten, yaitu Kabupaten Kepulauan Meranti yang mengalami perkembangan di tahun 2020 dari kluster dengan tingkat kemiskinan tertinggi menjadi kluster dengan tingkat kemiskinan sedang dan tidak pernah mengalami penurunan maupun perkembangan lagi setelah itu.
Selanjutnya, kota/kabupaten yang dari sebelum masa COVID-19 sampai masa akhir COVID-19 selalu di kluster dengan tingkat kemiskinan tertinggi terdapat pada Tabel berikut.
# Selalu di golongan termiskin
cluster_2019_to_2022%>% filter(c2022==3&c2019==3&c2020==3&c2021==3)## KOTKAB c2019 c2020 c2021 c2022
## 1 Lombok Utara 3 3 3 3
## 2 Sumba Timur 3 3 3 3
## 3 Timor Tengah Selatan 3 3 3 3
## 4 Sumba Tengah 3 3 3 3
## 5 Sabu Raijua 3 3 3 3
## 6 Maluku Barat Daya 3 3 3 3
## 7 Teluk Wondama 3 3 3 3
## 8 Teluk Bintuni 3 3 3 3
## 9 Tambrauw 3 3 3 3
## 10 Maybrat 3 3 3 3
## 11 Manokwari Selatan 3 3 3 3
## 12 Pegunungan Arfak 3 3 3 3
## 13 Jayawijaya 3 3 3 3
## 14 Nabire 3 3 3 3
## 15 Paniai 3 3 3 3
## 16 Puncak Jaya 3 3 3 3
## 17 Mappi 3 3 3 3
## 18 Asmat 3 3 3 3
## 19 Yahukimo 3 3 3 3
## 20 Pegunungan Bintang 3 3 3 3
## 21 Tolikara 3 3 3 3
## 22 Waropen 3 3 3 3
## 23 Supiori 3 3 3 3
## 24 Mamberamo Raya 3 3 3 3
## 25 Nduga 3 3 3 3
## 26 Lanny Jaya 3 3 3 3
## 27 Mamberamo Tengah 3 3 3 3
## 28 Yalimo 3 3 3 3
## 29 Puncak 3 3 3 3
## 30 Dogiyai 3 3 3 3
## 31 Intan Jaya 3 3 3 3
## 32 Deiyai 3 3 3 3Berdasarkan Tabel di atas, ada 32 kota/kabupaten yang selalu berada di kluster dengan tingkat kemiskinan tertinggi. Berdasarkan fakta bahwa 32 kota/kabupaten ini selalu berada di kluster dengan tingkat kemiskinan tertinggi, maka 32 kota/kabupaten ini merupakan kota/kabupaten yang membutuhkan perhatian lebih dari pemerintah.
Data Mapping
Adapun pemetaan untuk kota/kabupaten yang selalu berada di kluster dengan tingkat kemiskinan tertinggi dapat dilihat pada Gambar berikut.
shp <- st_read("/vsicurl/https://raw.githubusercontent.com/juliansalomo/dataset/main/BATAS%20KABUPATEN%20KOTA%20DESEMBER%202019%20DUKCAPIL/BATAS%20KABUPATEN%20KOTA%20DESEMBER%202019%20DUKCAPIL.shp")## Reading layer `BATAS%20KABUPATEN%20KOTA%20DESEMBER%202019%20DUKCAPIL' from data source `/vsicurl/https://raw.githubusercontent.com/juliansalomo/dataset/main/BATAS%20KABUPATEN%20KOTA%20DESEMBER%202019%20DUKCAPIL/BATAS%20KABUPATEN%20KOTA%20DESEMBER%202019%20DUKCAPIL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 515 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -11.00762 xmax: 141.0194 ymax: 6.07694
## Geodetic CRS: WGS 84DICT <- read.csv2("https://raw.githubusercontent.com/juliansalomo/dataset/main/dict.csv")
shp$id <- (0:514)%>%as.numeric() %>% as.character()
shpnew <- shp %>% filter(id!=0)
shpnew$KOTKAB <- with(DICT, KOTKAB[match(shpnew$id, id)])#Mapping the Persistence of Highest Poverty Clusters in Cities/Districts from 2019 to 2022
filt1 <- cluster_2019_to_2022 %>% filter(c2019*c2020*c2021*c2022==3^4)
KLASTER3 <- inner_join(shpnew, dataprofile2019[,c(1,13)],by = "KOTKAB")%>% filter(KOTKAB %in% filt1$KOTKAB)
# Load your sf dataset
sf_data <- KLASTER3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping the Persistence of Highest Poverty Clusters in Cities/Districts from 2019 to 2022") +
theme_minimal()
# Show the map
map
Mapping the Persistence of Highest Poverty Clusters in Cities/Districts from 2019 to 2022
Dari Gambar di atas bisa dikatakan mayoritas kota/kabupaten dengan tingkat kemiskinan tertinggi paling banyak berada di Pulau Papua, dan sisanya tersebar di pulau lain.
Selanjutnya kota/kabupaten yang masuk ke kluster dengan tingkat kemiskinan tertinggi pada akhir periode karena efek COVID-19 adalah sebagai berikut
# Mapping Clusters of Cities/Districts Classified as Poorest after COVID-19
filt2 <- cluster_2019_to_2022 %>% filter(c2019!=3 & c2022==3)
filt2## KOTKAB c2019 c2020 c2021 c2022
## 1 Sumba Barat 2 2 3 3
## 2 Lembata 2 2 3 3
## 3 Rote Ndao 2 2 3 3
## 4 Manggarai Timur 2 2 3 3
## 5 Maluku Tenggara Barat 2 3 3 3KLASTER3 <- inner_join(shpnew, dataprofile2022[,c(1,13)],by = "KOTKAB")%>% filter(KOTKAB %in% filt2$KOTKAB)
# Load your sf dataset
sf_data <- KLASTER3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping Clusters of Cities/Districts Classified as Poorest after COVID-19") +
theme_minimal()
# Show the map
map
Mapping Clusters of Cities/Districts Classified as Poorest after COVID-19
Selanjutnya pemetaan kota/kabupaten yang mengalami peningkatan tingkat kemiskinan sesudah COVID-19 dapat dilihat pada Gambar berikut.
# Mapping Clusters of Cities/Districts with Increasing Poverty Rates from 2019 to 2022
filt0 <- cluster_2019_to_2022%>% filter(c2019<c2022)
KLASTER3 <- inner_join(shpnew, dataprofile2022[,c(1,13)],by = "KOTKAB")%>% filter(KOTKAB %in% filt0$KOTKAB)
# Load your sf dataset
sf_data <- KLASTER3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = c("#95cb05","#ed0302"),name = "Cluster") +
labs(title = "Mapping Clusters of Cities/Districts with Increasing Poverty Rates from 2019 to 2022") +
theme_minimal()
# Show the map
map
Mapping Clusters of Cities/Districts with Increasing Poverty Rates from 2019 to 2022
Selanjutnya kota/kabupaten yang masuk ke kluster dengan tingkat kemiskinan tertinggi di tahun 2019 terdapat pada Tabel berikut, dan pemetaannya dapat dilihat pada Gambar setelahnya.
#Mapping Clusters with the Highest Poverty Rates in 2019
filt3 <- dataprofile2019[,c(1,13)] %>% filter(cluster == 3)
filt3## KOTKAB cluster
## 1 Kepulauan Meranti 3
## 2 Lombok Utara 3
## 3 Sumba Timur 3
## 4 Timor Tengah Selatan 3
## 5 Sumba Tengah 3
## 6 Sabu Raijua 3
## 7 Maluku Barat Daya 3
## 8 Teluk Wondama 3
## 9 Teluk Bintuni 3
## 10 Sorong 3
## 11 Tambrauw 3
## 12 Maybrat 3
## 13 Manokwari Selatan 3
## 14 Pegunungan Arfak 3
## 15 Jayawijaya 3
## 16 Nabire 3
## 17 Kepulauan Yapen 3
## 18 Paniai 3
## 19 Puncak Jaya 3
## 20 Boven Digoel 3
## 21 Mappi 3
## 22 Asmat 3
## 23 Yahukimo 3
## 24 Pegunungan Bintang 3
## 25 Tolikara 3
## 26 Waropen 3
## 27 Supiori 3
## 28 Mamberamo Raya 3
## 29 Nduga 3
## 30 Lanny Jaya 3
## 31 Mamberamo Tengah 3
## 32 Yalimo 3
## 33 Puncak 3
## 34 Dogiyai 3
## 35 Intan Jaya 3
## 36 Deiyai 3idn_cluster2019_3 <- inner_join(shpnew, filt3,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2019_3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping Clusters with the Highest Poverty Rates in 2019") +
theme_minimal()
# Show the map
map
Mapping Clusters with the Highest Poverty Rates in 2019
Selanjutnya kota/kabupaten yang masuk ke kluster dengan tingkat kemiskinan tertinggi di tahun 2020 terdapat pada Tabel berikut, dan pemetaannya dapat dilihat pada Gambar setelahnya.
#Mapping Clusters with the Highest Poverty Rates in 2020
filt4 <- dataprofile2020[,c(1,13)] %>% filter(cluster == 3)
filt4## KOTKAB cluster
## 1 Lombok Utara 3
## 2 Sumba Timur 3
## 3 Timor Tengah Selatan 3
## 4 Sumba Tengah 3
## 5 Sabu Raijua 3
## 6 Maluku Tenggara Barat 3
## 7 Maluku Barat Daya 3
## 8 Teluk Wondama 3
## 9 Teluk Bintuni 3
## 10 Tambrauw 3
## 11 Maybrat 3
## 12 Manokwari Selatan 3
## 13 Pegunungan Arfak 3
## 14 Jayawijaya 3
## 15 Nabire 3
## 16 Kepulauan Yapen 3
## 17 Paniai 3
## 18 Puncak Jaya 3
## 19 Mappi 3
## 20 Asmat 3
## 21 Yahukimo 3
## 22 Pegunungan Bintang 3
## 23 Tolikara 3
## 24 Waropen 3
## 25 Supiori 3
## 26 Mamberamo Raya 3
## 27 Nduga 3
## 28 Lanny Jaya 3
## 29 Mamberamo Tengah 3
## 30 Yalimo 3
## 31 Puncak 3
## 32 Dogiyai 3
## 33 Intan Jaya 3
## 34 Deiyai 3idn_cluster2020_3 <- inner_join(shpnew, filt4,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2020_3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping Clusters with the Highest Poverty Rates in 2020") +
theme_minimal()
# Show the map
map
Mapping Clusters with the Highest Poverty Rates in 2020
Selanjutnya kota/kabupaten yang masuk ke kluster dengan tingkat kemiskinan tertinggi di tahun 2021 terdapat pada Tabel berikut, dan pemetaannya dapat dilihat pada Gambar setelahnya.
#Mapping Clusters with the Highest Poverty Rates in 2021
filt5 <- dataprofile2021[,c(1,13)] %>% filter(cluster == 3)
filt5## KOTKAB cluster
## 1 Lombok Utara 3
## 2 Sumba Barat 3
## 3 Sumba Timur 3
## 4 Timor Tengah Selatan 3
## 5 Lembata 3
## 6 Rote Ndao 3
## 7 Sumba Tengah 3
## 8 Manggarai Timur 3
## 9 Sabu Raijua 3
## 10 Maluku Tenggara Barat 3
## 11 Kepulauan Aru 3
## 12 Maluku Barat Daya 3
## 13 Teluk Wondama 3
## 14 Teluk Bintuni 3
## 15 Sorong 3
## 16 Tambrauw 3
## 17 Maybrat 3
## 18 Manokwari Selatan 3
## 19 Pegunungan Arfak 3
## 20 Jayawijaya 3
## 21 Nabire 3
## 22 Paniai 3
## 23 Puncak Jaya 3
## 24 Boven Digoel 3
## 25 Mappi 3
## 26 Asmat 3
## 27 Yahukimo 3
## 28 Pegunungan Bintang 3
## 29 Tolikara 3
## 30 Waropen 3
## 31 Supiori 3
## 32 Mamberamo Raya 3
## 33 Nduga 3
## 34 Lanny Jaya 3
## 35 Mamberamo Tengah 3
## 36 Yalimo 3
## 37 Puncak 3
## 38 Dogiyai 3
## 39 Intan Jaya 3
## 40 Deiyai 3idn_cluster2021_3 <- inner_join(shpnew, filt5,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2021_3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping Clusters with the Highest Poverty Rates in 2021") +
theme_minimal()
# Show the map
map
Mapping Clusters with the Highest Poverty Rates in 2021
Selanjutnya kota/kabupaten yang masuk ke kluster dengan tingkat kemiskinan tertinggi di tahun 2022 terdapat pada Tabel berikut, dan pemetaannya dapat dilihat pada Gambar setelahnya.
#Mapping Clusters with the Highest Poverty Rates in 2022
filt6 <- dataprofile2022[,c(1,13)] %>% filter(cluster == 3)
filt6## KOTKAB cluster
## 1 Lombok Utara 3
## 2 Sumba Barat 3
## 3 Sumba Timur 3
## 4 Timor Tengah Selatan 3
## 5 Lembata 3
## 6 Rote Ndao 3
## 7 Sumba Tengah 3
## 8 Manggarai Timur 3
## 9 Sabu Raijua 3
## 10 Maluku Tenggara Barat 3
## 11 Maluku Barat Daya 3
## 12 Teluk Wondama 3
## 13 Teluk Bintuni 3
## 14 Sorong 3
## 15 Tambrauw 3
## 16 Maybrat 3
## 17 Manokwari Selatan 3
## 18 Pegunungan Arfak 3
## 19 Jayawijaya 3
## 20 Nabire 3
## 21 Kepulauan Yapen 3
## 22 Paniai 3
## 23 Puncak Jaya 3
## 24 Boven Digoel 3
## 25 Mappi 3
## 26 Asmat 3
## 27 Yahukimo 3
## 28 Pegunungan Bintang 3
## 29 Tolikara 3
## 30 Waropen 3
## 31 Supiori 3
## 32 Mamberamo Raya 3
## 33 Nduga 3
## 34 Lanny Jaya 3
## 35 Mamberamo Tengah 3
## 36 Yalimo 3
## 37 Puncak 3
## 38 Dogiyai 3
## 39 Intan Jaya 3
## 40 Deiyai 3idn_cluster2022_3 <- inner_join(shpnew, filt6,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2022_3
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = shpnew, fill = "#487a47", color = NA) +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = "#ed0302",name = "Cluster") +
labs(title = "Mapping Clusters with the Highest Poverty Rates in 2022") +
theme_minimal()
# Show the map
map
Mapping Clusters with the Highest Poverty Rates in 2022
Adapun pemetaan untuk hasil kluster kota/kabupaten dari tahun 2019 sampai 2022 dapat dilihat pada gambar-gambar berikut
# Mapping poverty clustering in 2019
idn_cluster2019 <- inner_join(shpnew, dataprofile2019,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2019
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = c("#05b67f","#95cb05","#ed0302"),name = "Cluster") +
labs(title = "Mapping poverty clustering in 2019") +
theme_minimal()
# Show the map
map
# Mapping poverty clustering in 2020
idn_cluster2020 <- inner_join(shpnew, dataprofile2020,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2020
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = c("#05b67f","#95cb05","#ed0302"),name = "Cluster") +
labs(title = "Mapping poverty clustering in 2020") +
theme_minimal()
# Show the map
map
# Mapping poverty clustering in 2021
idn_cluster2021 <- inner_join(shpnew, dataprofile2021,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2021
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = c("#05b67f","#95cb05","#ed0302"),name = "Cluster") +
labs(title = "Mapping poverty clustering in 2021") +
theme_minimal()
# Show the map
map
# Mapping poverty clustering in 2022
idn_cluster2022 <- inner_join(shpnew, dataprofile2022,by = "KOTKAB")
# Load your sf dataset
sf_data <- idn_cluster2022
# Plot the map by cluster
map <- ggplot() +
geom_sf(data = sf_data, aes(fill = as.factor(cluster)), color = NA) +
scale_fill_manual(values = c("#05b67f","#95cb05","#ed0302"),name = "Cluster") +
labs(title = "Mapping poverty clustering in 2022") +
theme_minimal()
# Show the map
map