Kerentanan Sosial Terhadap COVID-19 di DKI Jakarta Menggunakan Analisis Faktor dan Kluster

Multivariate Analysis

Multivariate analysis adalah sekelompok metode statistik yang berfokus pada mempelajari banyak variabel bersama-sama sambil berfokus pada variasi yang dimiliki oleh variabel-variabel tersebut.

Multivariate Statistics deals with the treatment of data sets with a large number of dimensions.

Analisis Faktor / Factor Analysis

Analisis faktor bertujuan untuk menerangkan struktur ragam-peragam melalui kombinasi linier dari variabel-variabel pembentuknya. Pada dasarnya, analisis faktor mendekatkan data pada suatu pengelompokkan berdasarkan adanya keeratan hubungan antar-dimensi pembentuk faktor. Statistik kunci yang relevan dengan analisis faktor antara lain

Principal Component Analysis (PCA)

PCA merupakan suatu algoritma unsupervised learning yang bertujuan untuk dimentionality reduction. Salah satu tujuan dari PCA adalah untuk mengubah variabel asli yang saling berkorelasi satu sama lainnya menjadi satu set variabel baru yang memiliki korelasi lebih kecil atau saling bebas.

Ide dasar dari PCA adalah untuk membuat sumbu (axis) baru yang dapat menangkap informasi sebesar mungkin. Sumbu baru ini adalah yang dinamakan sebagai Principal Component (PC). Untuk melakukan dimensionality reduction, kita akan memilih beberapa PC untuk dapat merangkum informasi yang dibutuhkan. Nilai masing-masing PC dibentuk dari matriks varian-kovarians dari suatu set variabel melalui kombinasi linear dari variabel-variabel tersebut.

\[ PC1 = a_{11}x_1 + a_{12}x_2 + ... + a_{1p}x_p \] di mana:

• vektor \(a_1\) adalah eigen vector dari suatu matriks varian kovarians

• Nilai \(x_i\) adalah nilai suatu data pada baris ke-i

Hal-hal yang menjadi pedoman dalam menentukan jumlah PC, yaitu:

• 1. Mempertahankan PC yang cukup untuk menjelaskan persentase yang spesifik dari total variansi/informasi
• 2. Mempertahankan PC yang memiliki eigen value lebih besar dibandingkan rata-rata eigen value dari setiap PC

Fuzzy C-Means Clustering

Algoritma Fuzzy C-Means (FCM) pertama kali diperkenalkan oleh Dunn pada tahun 1973 dan dikembangkan kembali oleh Bezdek pada tahun 1981 sebagai metode clustering yang menggunakan konsep fuzzy.

Clustering adalah pengelompokan data menjadi beberapa cluster, di mana setiap data pada cluster yang sama akan memiliki kemiripan dibandingkan dengan cluster yang berbeda.

FCM adalah algoritma partitial clustering sehingga langkah pertama adalah menentukan banyaknya jumlah cluster pada data. Dilanjutkan dengan mencari pusat cluster (centroid), yang akan menandai posisi rata-rata untuk masing-masing cluster, lalu diperbarui setiap iterasi dan terakhir data dikelompokkan berdasarkan jarak antar data dengan setiap centroid.

Uji Kruskal Wallis

Pengujian ini merupakan uji non parametrik yang digunakan untuk menentukan ada atau tidaknya perbedaan yang signifikan pada dua atau lebih kelompok. Asumsi yang diperlukan dalam melakukan pengujian Kruskal Wallis adalah:

• Variabel dependen berskala numerik baik interval/rasio atau skala ordinal
• Variabel independent berskala kategorik dengan 2 atau lebih kategorik
• Independen. Artinya sampel di tiap kategori harus bebas satu sama lain.

Melakukan pengecekan missing value

melkukan pengecekan missing value pada data untuk setiap variabel.

colSums(is.na(rentan))

#>    Kecamatan    perempuan KK_perempuan       lansia pengangguran    wajib_ktp 
#>            0            0            0            0            0            0 
#>       banjir     informal           RS    puskesmas  positiv_cov        rawat 
#>            0            0            0            0            0            0 
#>       isoman         Kode 
#>            0            0

1. Gambaran umum variabel penyusun indeks kerentanan sosial

Ini adalah langkah yang baik untuk memahami data terlebih dahulu dan mencoba mengumpulkan sebanyak mungkin wawasan darinya.

Exploratory Data Analysis refers to the critical process of performing initial investigations on data so as to discover patterns,to spot anomalies,and to check the summary statistics and graphical representations

summary(rentan)

#>   Kecamatan           perempuan      KK_perempuan       lansia      
#>  Length:44          Min.   :49.11   Min.   :14.78   Min.   : 5.873  
#>  Class :character   1st Qu.:49.43   1st Qu.:18.93   1st Qu.: 7.844  
#>  Mode  :character   Median :49.67   Median :21.11   Median : 8.845  
#>                     Mean   :49.73   Mean   :21.09   Mean   : 8.997  
#>                     3rd Qu.:49.98   3rd Qu.:22.94   3rd Qu.: 9.776  
#>                     Max.   :50.98   Max.   :27.49   Max.   :13.219  
#>   pengangguran      wajib_ktp         banjir         informal     
#>  Min.   :0.1950   Min.   :67.78   Min.   : 0.00   Min.   : 6.148  
#>  1st Qu.:0.3217   1st Qu.:71.57   1st Qu.: 0.00   1st Qu.: 9.543  
#>  Median :0.5162   Median :72.66   Median : 9.00   Median :10.540  
#>  Mean   :0.5611   Mean   :72.78   Mean   :12.86   Mean   :12.412  
#>  3rd Qu.:0.6276   3rd Qu.:73.73   3rd Qu.:18.75   3rd Qu.:16.318  
#>  Max.   :1.7709   Max.   :79.05   Max.   :59.00   Max.   :25.738  
#>        RS          puskesmas      positiv_cov        rawat       
#>  Min.   :0.000   Min.   : 3.00   Min.   :  242   Min.   :  0.00  
#>  1st Qu.:1.000   1st Qu.: 6.00   1st Qu.: 5698   1st Qu.: 49.75  
#>  Median :1.500   Median : 7.50   Median : 7966   Median : 62.50  
#>  Mean   :1.818   Mean   : 7.75   Mean   : 7587   Mean   : 65.61  
#>  3rd Qu.:3.000   3rd Qu.: 9.25   3rd Qu.: 9622   3rd Qu.: 77.00  
#>  Max.   :7.000   Max.   :14.00   Max.   :14329   Max.   :142.00  
#>      isoman           Kode          
#>  Min.   :  0.00   Length:44         
#>  1st Qu.: 41.75   Class :character  
#>  Median : 62.50   Mode  :character  
#>  Mean   : 67.93                     
#>  3rd Qu.: 88.50                     
#>  Max.   :175.00

Berdasarkan hasil summary(rentan), secara umum terdapat 49,73 persen jumlah penduduk perempuan yang ditunjukan oleh rata-rata perempuan dengan 21,09 persen merupakan kepala keluarga yang berjenis kelamin perempuan dari variabel KK_perempuan. Sedangkan untuk rata-rata jumlah penduduk lansia di Jakarta sebesar 8,997 persen yang ditunjukan oleh lansia.

rmarkdown::paged_table(rentan[which(rentan$pengangguran < 0.1951 | rentan$pengangguran > 1.770), ])

Kemudian untuk variabel pengangguran, secara rata-rata sebesar 0,5611 persen dengan nilai maksimum 1,77 persen berada di Kecamatan Kep. Seribu Selatan dan nilai terendah 0,195 berada di Kecamatan Kembangan. Perbedaan antara nilai maksimum dan minimum yang hampir 10 kali lipat ini menunjukkan belum meratanya tingkat pengangguran di DKI Jakarta. Selain itu, terdapat 12,412 persen pekerja informal di DKI Jakarta.

Membahas lebih lanjut lagi mengenai variabel yang digunakan dalam penelitian ini, KTP merupakan salah satu dokumen administrasi yang sangat penting dan harus dimiliki seluruh masyarakat Indonesia. Tanpa adanya kepemilikan dokumen ini, maka akan terjadi kesulitan dalam menggunakan berbagai layanan publik, seperti kesehatan, perbankan, pendidikan, pendataan pemerintah, dan lain sebagainya. Namun, persentase penduduk wajib KTP di DKI Jakarta yang ditunjukkan oleh wajib_KTP secara umum hanya sebesar 72,78 persen, yang menandakan belum semua masyarakat sadar akan pentingnya kepemilikan KTP.

Dari hasil summary(), diketahui bahwa setiap variabel memiliki range nilai yang cukup berbeda, hal ini menyebabkan data yang dimiliki perlu untuk dilakukan scaling/standardisasi.

2. Melihat visualisasi antar beberapa variabel

Pertama dilakukan pemeriksaan nilai korelasi antar variabel dengan melihat scatter plot antar masing-masing nilai variabel berpasangan pada data.

ggpairs(rentan[,c(2:13)], showStrips = F, upper = list(continuous = wrap("cor", size=2))) + 
  theme(axis.text = element_text(colour = "black", size = 7),
        strip.background = element_rect(fill = "#6e8f71"),
        strip.text = element_text(colour = "white", size = 6,
                                  face = "bold"))

Kemudian dengan menggunakan bantuan ggplot ingin dicari informasi mengenai persentase lansia pada masing-masing kecamatan di provinsi DKI Jakarta.

ggplot(data = rentan, mapping = aes(x=  lansia, y= reorder(Kecamatan,lansia))) +
  geom_col(aes(fill = lansia))  + # untuk membuat barplot 
  labs(
    title = "Persentase Jumlah Lansia pada masing-masing Kecamatan di Provinsi DKI Jakarta",
    x = "Persentase Lansia",
    y = "Kecamatan"
  ) +
  scale_fill_gradient(low = "#a9d1b1", high = "#4c614e") +
  geom_text(mapping = aes(label=round(lansia, digits = 4)), # menambahkan informasi label masing-masing kategori dengan variasi box
            col = "white", # memberikan warna pada text
            nudge_x = -1, # menggeser text berdasarkan sumbu x
            label.size = 10) + 
  theme_minimal() +
  theme(legend.position = "none") # untuk menghilangkan legend.

Dari plot di atas dapat ditarik insight bahwa Kecamatan Gambir merupakan kecamatan dengan persentase lansia terbanyak yaitu sebesar 13.22%. Diikuti dengan Kecamatan Senen sebesar 12.16%. Kecamatand dengan persentase lansia terkecil terdapat pada Kecamatan Kalideres.

Melakukan uji multikolinearitas pada seluruh variabel.

Tujuan dari membentuk plot korelasi adalah untuk menguji kolinearitas dari data. Seluruh variabel harus memiliki korelasi, tetapi tidak terlalu tinggi.

ggcorr(rentan, label = T, hjust = 1, layout.exp = 2, size = 3, label_size = 2, label_round = 2, low = "#d7fc5b", mid = "white", high = "#5d8064")

Correlation plot di atas yang digunakan untuk menguji multikolinearitas dari data. Pada project ini, digunakan batas (>0,85) sebagai batas korelasi ekstrem terhadap setiap variabel. Dapat dilihat bahwa semua variabel tidak memiliki nilai korelasi yang lebih dari 0,85.

Kaiser–Meyer–Olkin (KMO) - Measure of Sampling Adequacy(MSA)

Untuk dapat melakukan pengujian KMO dan MSA, pada R disediakan suatu function KMO() dari library psych.

rentan_cor <- cor(rentan[,2:13])
rentan_KMO<-KMO(rentan_cor)

Melakukan inspeksi nilai KMO

rentan_KMO$MSA

#> [1] 0.7315378

Untuk nilai KMO didapatkan sebesar 0,7315378. Syarat dapat dilanjutkannya proses analisis faktor adalah nilai KMO lebih besar dari nilai 0,5. Dari hasil tersebut dapat dikatakan bahwa proses analisis faktor bisa dilanjutkan ke proses selanjutnya karena syarat dilanjutkannya proses analisis faktor terpenuhi.

Melakukan inspeksi nilai MSA masing-masing variabel

rmarkdown::paged_table(as.data.frame(rentan_KMO$MSAi))

Dari hasil output nilai MSA masing-masing variabel di atas didapatkan seluruh nilai MSA dari setiap variabel (> 0,5). Sehingga bisa dinyatakan bahwa semua variabel dapat diprediksi dan dianalisis lebih lanjut

The Bartlett’s Test statistics

Bartlett’s test digunakan untuk mengkonfirmasi bahwa data berkorelasi secara berpasangan (Mickey & Sharma, 1997). Dengan kata lain uji bartlett digunakan untuk mengetahui homogenitas variance.

• Hipotesis uji
  • H0 = Varians di antara masing-masing variabel adalah sama (variance homogen)
  • H1 = Setidaknya satu variabel memiliki varians yang tidak sama dengan yang lain (variance tidak homogen)
• Pengujian yang diharapkan adalah data memiliki variance yang tidak homogen (tolak H0).

bartlett.test(rentan[,2:13])

#> 
#>  Bartlett test of homogeneity of variances
#> 
#> data:  rentan[, 2:13]
#> Bartlett's K-squared = 5144, df = 11, p-value < 0.00000000000000022

Hasil perhitungan dengan R dihasilkan nilai Bartlett test of homogeneity variances sebesar 5144 dengan signifikansi sebesar 0,000. Dengan demikian, Bartlett test of homogeneity of variances memenuhi persyaratan karena signifikansi di bawah 0,05 (5%), sehingga variance tidak homogen.

Membentuk faktor-faktor penyusun kerentanan sosial

Untuk melakukan exploratory factor analysis menggunakan MinRes (Minimum Residual), fungsi yang tersedia di R adalah fa() dari library psych.

Parameter fungsi fa() - r: Correlation/covariance matrix dari data awal

fa_rentan <- fa(r = rentan_cor)
fa_rentan

#> Factor Analysis using method =  minres
#> Call: fa(r = rentan_cor)
#> Standardized loadings (pattern matrix) based upon correlation matrix
#>                MR1    h2   u2 com
#> perempuan    -0.11 0.013 0.99   1
#> KK_perempuan -0.59 0.352 0.65   1
#> lansia       -0.68 0.459 0.54   1
#> pengangguran -0.16 0.026 0.97   1
#> wajib_ktp    -0.58 0.334 0.67   1
#> banjir        0.51 0.260 0.74   1
#> informal     -0.69 0.480 0.52   1
#> RS            0.13 0.018 0.98   1
#> puskesmas     0.63 0.402 0.60   1
#> positiv_cov   0.85 0.721 0.28   1
#> rawat         0.64 0.406 0.59   1
#> isoman        0.72 0.519 0.48   1
#> 
#>                 MR1
#> SS loadings    3.99
#> Proportion Var 0.33
#> 
#> Mean item complexity =  1
#> Test of the hypothesis that 1 factor is sufficient.
#> 
#> The degrees of freedom for the null model are  66  and the objective function was  9.27
#> The degrees of freedom for the model are 54  and the objective function was  5.98 
#> 
#> The root mean square of the residuals (RMSR) is  0.23 
#> The df corrected root mean square of the residuals is  0.26 
#> 
#> Fit based upon off diagonal values = 0.66
#> Measures of factor score adequacy             
#>                                                    MR1
#> Correlation of (regression) scores with factors   0.96
#> Multiple R square of scores with factors          0.91
#> Minimum correlation of possible factor scores     0.83

Secara default, objek factor analysis yang terbentuk menghasilkan sejumlah 1 faktor. Jumlah ini dapat saja belum optimal, sehingga perlu untuk mencari jumlah faktor yang optimal dari data yang kita miliki.

Salah satu cara untuk mencari jumlah faktor yang optimal dari data, digunakan sebuah fungsi yaitu fa.parallel(). Metode ini mencari jumlah faktor optimal dari nilai eigen yang berurutan, kemudian penurunan tajam dalam plot menunjukkan jumlah komponen atau faktor yang sesuai untuk diekstraksi.

nfactors <- fa.parallel(rentan[,2:13], plot = T)

#> Parallel analysis suggests that the number of factors =  2  and the number of components =  2

summary(nfactors)

Dari hasil paralell analysis scree plot di atas dihasilkan factors=2 sebagai jumlah faktor yang optimal pada data rentan.

Selanjutnya akan kembali dilakukan analisis faktor menggunakan jumlah faktor optimal yang sudah didapatkan. Parameter tambahan fungsi fa() - nfactors: banyaknya faktor yang akan dibentuk - rotate : tipe rotasi. Rotasi meminimalkan kompleksitas faktor untuk membuat struktur lebih sederhana untuk ditafsirkan.

# jumlah faktor optimal yang didapat
nfact <- nfactors$nfact

#factor analysis
fa_rentan_opt <- fa(r = rentan_cor, nfactors = nfact, rotate = "varimax")
fa_rentan_opt

#> Factor Analysis using method =  minres
#> Call: fa(r = rentan_cor, nfactors = nfact, rotate = "varimax")
#> Standardized loadings (pattern matrix) based upon correlation matrix
#>                MR1   MR2   h2      u2 com
#> perempuan     0.22  0.48 0.27 0.72599 1.4
#> KK_perempuan -0.15  0.83 0.71 0.29360 1.1
#> lansia       -0.21  0.91 0.88 0.12095 1.1
#> pengangguran -0.50 -0.37 0.38 0.61504 1.8
#> wajib_ktp    -0.08  0.93 0.87 0.12878 1.0
#> banjir        0.50 -0.17 0.28 0.72403 1.2
#> informal     -0.36  0.66 0.56 0.43507 1.5
#> RS            0.38  0.27 0.22 0.78159 1.8
#> puskesmas     0.61 -0.22 0.42 0.57622 1.3
#> positiv_cov   0.99 -0.17 1.00 0.00042 1.1
#> rawat         0.88  0.03 0.77 0.22694 1.0
#> isoman        0.80 -0.15 0.66 0.33818 1.1
#> 
#>                        MR1  MR2
#> SS loadings           3.65 3.39
#> Proportion Var        0.30 0.28
#> Cumulative Var        0.30 0.59
#> Proportion Explained  0.52 0.48
#> Cumulative Proportion 0.52 1.00
#> 
#> Mean item complexity =  1.3
#> Test of the hypothesis that 2 factors are sufficient.
#> 
#> The degrees of freedom for the null model are  66  and the objective function was  9.27
#> The degrees of freedom for the model are 43  and the objective function was  1.6 
#> 
#> The root mean square of the residuals (RMSR) is  0.05 
#> The df corrected root mean square of the residuals is  0.07 
#> 
#> Fit based upon off diagonal values = 0.98

Dari output di atas dengan menggunakan faktor sebanyak 2 didapatkan beberapa informasi sebagai berikut:

• SS Loadings: Nilai Eigen.

• Proportion Var: Memberitahu kita berapa banyak varians/informasi keseluruhan yang diperhitungkan faktor dari semua variabel.

• Cumulative Var: Nilai akumulasi dari proportion var.

• Proportion Explained: Jumlah variansi relatif yang terjelaskan. Merupakan pembagian dari proportion var/ sum(proportion var).

• Cumulative Proportion: Jumlah akumulasi dari proportion explained.

• MR: Minimum Residual/factor dapat disebut juga sebagai komponen utama.

• h2: jumlah varians dalam item/variabel yang dijelaskan oleh faktor-faktor (yang dipertahankan). Ekuivalen dengan SS Loadings untuk variabel / communality

• u2: merupakan nilai 1-h2 adalah residual variance.

• com: item complexity.

Dengan menggunakan 2 faktor, maka informasi keseluruhan yang diperhitungkan faktor dari semua variabel sebanyak 0,59 atau 59%.

Faktor-Faktor Pembentuk Kerentanan Sosial di Provinsi DKI Jakarta

Salah satu tujuan dari factor analysis adalah melakukan ekstraksi dengan menghasilkan variabel laten. Variabel laten yang dihasilkan nantinya akan dijadikan sebagai faktor-faktor pembentuk kerentanan sosial pada data.

fa.diagram(fa_rentan_opt, rsize = 1)

Arah dari masing-masing variabel ditentukan berdasarkan kajian literatur dari (Rufat et al., 2015), dimana arah tersebut ditentukan berdasarkan pengaruhnya terhadap kerentanan sosial. Arah positif berpengaruh terhadap peningkatan kerentanan sosial, sebaliknya arah negatif menurunkan tingkat kerentanan sosial.

Dari plot di atas dilakukan pengekstrakan terhadap 2 faktor/ komponen utama, yaitu kesehatan dan pekerjaan serta kependudukan dan KRT perempuan.

• Faktor Kesehatan dan Pekerjaan

Faktor kesehatan terdiri dari tujuh variabel pembentuk, yaitu jumlah kasus positif COVID-19, jumlah penduduk yang dirawat akibat COVID-19, jumlah penduduk yang melakukan isolasi mandiri akibat COVID-19, jumlah puskesmas, persentase pengangguran, jumlah daerah rawan banjir, dan jumlah RS umum. Pada faktor ini menjelaskan 30 persen varians yang membentuk kerentanan sosial. Lima variabel berpengaruh positif dan dua variabel berpengaruh negatif dalam kerentanan sosial.

Daerah yang memiliki jumlah penderita COVID-19 yang tinggi meningkatkan kontribusi terhadap kerentanan, dimana penularan COVID-19 yang sangat cepat dapat menimbulkan daerah tersebut diberlakukan pembatasan kegiatan, sehingga memperlambat pemulihan terhadap dampak yang ditimbulkan COVID-19.

Kemudian, jumlah puskesmas yang banyak dapat menurunkan kerentanan sosial, dikarenakan dengan tersedianya fasilitas kesehatan yang memadai maka akan cepat melakukan penanganan dan pemulihan kesehatan kepada masyarakat yang terdampak sehingga pemulihan terhadap dampak COVID-19 dapat dipercepat.

• Faktor Kependudukan dan KRT Perempuan

Faktor kependudukan terdiri dari lima variabel pembentuk, yaitu persentase penduduk wajib KTP, persentase penduduk lansia, persentase kepala keluarga berjenis kelamin perempuan, persentase pekerja di sektor informal, dan persentase penduduk perempuan. Pada faktor ini menjelaskan 28 persen varians yang membentuk kerentanan sosial.

Struktur penduduk yang didominasi oleh penduduk lansia menyebabkan ketergantungan yang besar terhadap orang lain, sehingga ketika terjadi COVID-19 kelompok ini akan sangat rentan. Kelompok ini akan lebih sulit untuk memulihkan diri dikarenakan kemampuan fisik yang sudah menurun, sehingga mobilitas terhambat.

Selain itu, perempuan menjadi kelompok yang lebih sulit untuk pulih dari dampak COVID-19. Selanjutnya, masyarakat yang memiliki KTP akan lebih mudah mendapatkan akses terhadap pendataan, sehingga dengan pendataan yang baik dapat menyalurkan bantuan yang terarah dan terukur sesuai KTP dan menyebabkan pemulihan lebih cepat.

Indeks Kerentanan Sosial

Ada banyak metode untuk melakukan pembobotan dalam agregasi indeks. Peneliti seperti:

• (Cutter et al., 2003) menggunakan pembobot tertimbang,

• (Solangaarachchi et al., 2012) menggunakan pembobot berdasarkan kontribusi terhadap varians terjelaskan,

• (Siagian et al., 2014) menggunakan pembobot tidak sama, dan

• (Rygel et al., 2006) menggunakan pareto-ranking. Tidak ada kesepakatan yang membatasi penggunaan pembobotan yang berbeda.

Dalam project ini, digunakan pembobotan tidak sama untuk melakukan agregasi indeks. Pembobot didapatkan dari pembagian antara varians terjelaskan masing-masing komponen dibagi dengan total varians terjelaskan atau nilai proportion explained dari object fa_rentan_opt. Dengan cara ini komponen yang memiliki varians terjelaskan yang lebih tinggi berpengaruh lebih besar terhadap kerentanan sosial. Perhitungan kerentanan sosial mengikuti persamaan sebagai berikut:

\[SoVI=(0.52*komponen1)+(0.48*komponen2)\] SoVI merujuk peda Social Vulnarability Index.

Melakukan PCA pada data rentan

PCA bertujuan untuk melakukan reduksi data dari sekumpulan variabel (Mickey & Sharma, 1997). Pada analisis ini, PCA digunakan untuk menghitung nilai SoVI masing-masing kecamatan yang ada pada dataset. Untuk itu akan dilakukan inspeksi untuk pca$x dimana merupakan nilai hasil proyeksi titik ke PC untuk tiap baris. Digunakan untuk mendapatkan nilai data yang baru.

rentan_pca<-prcomp(rentan[,2:13], scale. = T)
rmarkdown::paged_table(as.data.frame(rentan_pca$x))

Sesuai dengan hasil dari factor analysis, akan kita ambil sebanyak 2 principal component untuk kebutuhan pemetaan SoVI masing-masing kecamatan di Provinsi DKI Jakarta.

pc_keep <- rentan_pca$x[, 1:2] %>% 
  as.data.frame()
rmarkdown::paged_table(pc_keep)

Melakukan perhitungan persamaan indeks kerentanan sosial.

SoVI <- pc_keep %>% 
  mutate(SoVI = 0.52*abs(PC1) + 0.48*abs(PC2))
rmarkdown::paged_table(SoVI)

Selanjutnya, dilakukan pengabungan nilai hasil perhitungan Social Vulnerability Index dengan nama-nama kecamatan dan kode area yang ada pada dataset rentan.

rentan_new <- rentan %>% 
  select_if(~!is.numeric(.)) %>% # ambil kolom Kecamatan
  cbind(SoVI = SoVI[,3]) # gabungkan dengan kolom SoVI
  
rmarkdown::paged_table(rentan_new)

Dari indeks kerentanan sosial yang terbentuk, kami melakukan pengkategorian indeks menjadi tiga, yaitu rendah, sedang, dan tinggi. Kami melakukan pengkategorian indeks menggunakan Natural Breaks pada R menggunakan function getJenksBreaks() dari library BAMMtools.

# Membagi SoVI ke dalam 3 kategori
natural_breaks <- function(df, var) {
  var_breaks <- BAMMtools::getJenksBreaks(df[[var]], k = 4)
  df[[paste0('Jenks_', var)]] <- findInterval(df[[var]], var_breaks)
  df
}
rentan_new <- natural_breaks(rentan_new, "SoVI")

# Membuat Kategori Rendah Sedang Tinggi
rentan_new$Kategori <- ordered(rentan_new$Jenks_SoVI,
levels = c(1,2,3),
labels = c("Rendah", "Sedang", "Tinggi"))
rmarkdown::paged_table(rentan_new)

# Melihat threshold untuk masing-masing kategori 
getJenksBreaks(rentan_new$SoVI, 4)

#> [1] 0.2734244 1.3269436 2.3094436 3.5813008

Hasil yang didapatkan yaitu 0.2734 sampai 1.3269 untuk kategori rendah; 1.3269 sampai 2.3094 untuk kategori sedang; dan 2.3094 sampai 3.5813 untuk kategori tinggi. Secara umum, sebagian besar masyarakat di wilayah DKI Jakarta memiliki kerentanan sosial terhadap COVID-19 yang berada pada rentang kategori sedang-rendah.

Visualisasi Geospasial Kerentanan Sosial Kecamatan di Provinsi DKI Jakarta

Pemetaan dalam bentuk visual dilakukan untuk melihat persebaran indeks kerentanan sosial menurut kecamatan di DKI Jakarta. Visualisasi yang akan dilakukan dinamakan visualisasi geospasial, kita bisa memberikan informasi yang hanya dengan bantuan peta maka informasi tersebut dapat tersampaikan dengan baik.

Pada project ini akan dibuat visualisasi geocoding berupa choropleth.

Pertama-tama dilakukan penarikan data SHP untuk keseluruhan provinsi di Indonesia yang diambil melalui website GADM. Kemudian melakukan subseting data dengan hanya mengambil Provinsi DKI Jakarta dengan bantuan MapShaper.

Langkah kedua, untuk membaca data .json di R, kita akan menggunakan package geojsonio. Pada fungsi geojson_read(), kita menentukan what = “sp” untuk memastikan data yang kita baca merupakan kelas spasial, bukan list. Kemudian untuk mengecek datanya terlebih dahulu, kita dapat fokus pada isi tabelnya dengan menggunakan @.

library(geojsonio)
jakarta_json <- geojson_read("dataset/gadm36_IDN_3.json", what = "sp")
rmarkdown::paged_table(jakarta_json@data)

Selanjutnya, kita melakukan convert object jakarta_json menjadi sebuah object sf. Object sf adalah sebuah dataframe yang berisi aspek spasial (geospasial).

jakarta_json_mod <- sf::st_as_sf(jakarta_json)

Langkah keempat setelah membuat object sf adalah melakukan data pre-processing untuk data tersebut. Beberapa diantaranya yaitu membuang kolom yang tidak berguna.

library(stringr)
jakarta_json_mod <- jakarta_json_mod %>% 
  # Agar mudah digabung
  mutate(NAME_3 = str_replace_all(NAME_3, fixed(" "), "") %>% str_to_title()) %>% 
  # melakukan pemotongan variabel/kolom yang tidak berguna untuk analisis
  dplyr::select(-c(id, NL_NAME_1, NL_NAME_2, NL_NAME_3, VARNAME_3, HASC_3, TYPE_3, ENGTYPE_3)) 

# melakukan inspeksi data dan tipenya
glimpse(jakarta_json_mod)

#> Rows: 47
#> Columns: 10
#> $ GID_0    <chr> "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN"…
#> $ NAME_0   <chr> "Indonesia", "Indonesia", "Indonesia", "Indonesia", "Indonesi…
#> $ GID_1    <chr> "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN.7…
#> $ NAME_1   <chr> "Jakarta Raya", "Jakarta Raya", "Jakarta Raya", "Jakarta Raya…
#> $ GID_2    <chr> "IDN.7.1_1", "IDN.7.1_1", "IDN.7.1_1", "IDN.7.1_1", "IDN.7.1_…
#> $ NAME_2   <chr> "Jakarta Barat", "Jakarta Barat", "Jakarta Barat", "Jakarta B…
#> $ GID_3    <chr> "IDN.7.1.1_1", "IDN.7.1.2_1", "IDN.7.1.3_1", "IDN.7.1.4_1", "…
#> $ NAME_3   <chr> "Cengkareng", "Grogolpetamburan", "Kalideres", "Kebonjeruk", …
#> $ CC_3     <chr> "3174070", "3174040", "3174080", "3174020", "3174010", "31740…
#> $ geometry <MULTIPOLYGON> MULTIPOLYGON (((106.7004 -6..., MULTIPOLYGON (((106.…

Langkah kelima adalah melakukan pengolahan data denganmencoba untuk menyocokkan isi dari data rentan yang nanti akan kita tampilkan informasinya di peta, dengan data spasial yang sebelumnya sudah kita baca (dan kita olah sedikit). Hal ini bermaksud agar tidak terjadinya bentrok dikarenakan nama daerah yang berbeda, dan agar lebih masuk akal dilakukan.

jakarta_rentan <- jakarta_json_mod %>% 
  left_join(rentan_new, by = c("GID_3" = "Kode")) %>% 
  arrange(NAME_3, NAME_0)
glimpse(jakarta_rentan)

#> Rows: 47
#> Columns: 14
#> $ GID_0      <chr> "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "IDN", "ID…
#> $ NAME_0     <chr> "Indonesia", "Indonesia", "Indonesia", "Indonesia", "Indone…
#> $ GID_1      <chr> "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN.7_1", "IDN…
#> $ NAME_1     <chr> "Jakarta Raya", "Jakarta Raya", "Jakarta Raya", "Jakarta Ra…
#> $ GID_2      <chr> "IDN.7.4_1", "IDN.7.2_1", "IDN.7.1_1", "IDN.7.3_1", "IDN.7.…
#> $ NAME_2     <chr> "Jakarta Timur", "Jakarta Pusat", "Jakarta Barat", "Jakarta…
#> $ GID_3      <chr> "IDN.7.4.1_1", "IDN.7.2.1_1", "IDN.7.1.1_1", "IDN.7.3.1_1",…
#> $ NAME_3     <chr> "Cakung", "Cempakaputih", "Cengkareng", "Cilandak", "Cilinc…
#> $ CC_3       <chr> "3172080", "3173050", "3174070", "3171030", "3175060", "317…
#> $ Kecamatan  <chr> "Cakung", "Cempaka Putih", "Cengkareng", "Cilandak", "Cilin…
#> $ SoVI       <dbl> 2.1771459, 1.5710675, 2.2455849, 0.6934557, 1.7166052, 1.57…
#> $ Jenks_SoVI <int> 2, 2, 2, 1, 2, 2, 1, NA, 3, 3, 2, 2, 1, 2, NA, 2, 2, 1, 1, …
#> $ Kategori   <ord> Sedang, Sedang, Sedang, Rendah, Sedang, Sedang, Rendah, NA,…
#> $ geometry   <MULTIPOLYGON> MULTIPOLYGON (((106.9508 -6..., MULTIPOLYGON (((10…

Akan tetapi yang perlu dinote bahwa, pada data spasial yang kita miliki terdapat 2 object yang redundan, di antaranya adalah kebayoran lama dengan variabel GID_3 = IDN.7.3.10_1 dan setiabudi dengan variabel GID_3 = IDN.7.3.3_1. Serta 1 object danausunterdll yang tidak terdapat pada data rentan. Sehingga nantinya ketiga object ini akan diabaikan/dikeluarkan.

jakarta_rentan_sf <- jakarta_rentan %>%  
  filter(GID_3 != "IDN.7.3.10_1") %>% 
  filter(GID_3 != "IDN.7.3.3_1") %>% 
  filter(GID_3 != "IDN.7.5.2_1") %>%
  sf::st_as_sf()

Dalam penyajian data, pertama-tama dimasukkan data rentan dalam leaflet. Lalu, akan ditentukan pewarnaan yang sesuai untuk choropleth yang nanti akan dibentuk. Palette ditentukan dari kuning hingga merah (silahkan lihat cheatsheet yang sudah disediakan). Selain itu, dibuat juga popup pada tiap tempat.

library(leaflet)
m <- leaflet(jakarta_rentan_sf)

popup.cont  <- paste("<h2><b>", jakarta_rentan_sf$NAME_3, "</b></h2>")
col <- jakarta_rentan_sf$SoVI
col <- col %>% replace(is.na(.), 0)
bins <- c(getJenksBreaks(jakarta_rentan_sf$SoVI, 4))
pal <- colorBin("YlOrRd", domain = col, bins = bins,)

Terakhir, kita akan membentuk peta dengan tiles yang disediakan oleh CartoDB menggunakan fungsi addProviderTiles. Lalu kita atur poligon dan legenda pada peta, serta mengubah nama tiap elemen pada legenda agar lebih masuk akal.

m %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(fillColor = pal(col),
              weight = 1,
              opacity = 1,
              color = "black",
              dashArray = "3",
              fillOpacity = 0.8,
              label = paste0("Kecamatan: ", jakarta_rentan_sf$NAME_3),
              popup = popup.cont) %>%
  addLegend("bottomright", 
            pal = pal,
            values = ~col,
            title = "Social Vulnerability Index",
            labFormat = labelFormat(digits = 2),
            opacity = 1)

Data pre-processing

Langkah pertama mendefinisikan data yang akan digunakan dalam analisis FCM ini ke dalam object rentan_FCMm dimana data ini berisi variabel-variabel bertipe data numerik.

rentan_FCM <- rentan[,2:13]
head(rentan_FCM)

#> # A tibble: 6 × 12
#>   perempuan KK_per…¹ lansia penga…² wajib…³ banjir infor…⁴    RS puske…⁵ posit…⁶
#>       <dbl>    <dbl>  <dbl>   <dbl>   <dbl>  <dbl>   <dbl> <dbl>   <dbl>   <dbl>
#> 1      49.2     18.1   6.22   0.532    70.2     31    6.50     1       9   11740
#> 2      50.1     23.0  10.9    0.538    75.1      0   13.9      2       4    4998
#> 3      49.3     18.5   6.59   0.516    70.8     55    6.79     1      10   13036
#> 4      50.4     20.8   9.31   0.267    73.7      2   10.9      3       6    8201
#> 5      49.5     20.5   5.98   0.653    68.9     38   11.3      1      10    7899
#> 6      49.6     17.1   7.00   0.871    71.2     15   14.3      0      11   10437
#> # … with 2 more variables: rawat <dbl>, isoman <dbl>, and abbreviated variable
#> #   names ¹​KK_perempuan, ²​pengangguran, ³​wajib_ktp, ⁴​informal, ⁵​puskesmas,
#> #   ⁶​positiv_cov

Melakukan assign nilai dari kolom Kecamatan yang ada pada data rentan kita menjadi nama baris.

rownames(rentan_FCM) <- rentan$Kecamatan

Normalisasi Min-Max

Teknik penskalaan Min-Max merupakan salah satu metode normaliasi dengan melakukan perubahan linear terhadap data asli. Teknik normalisasi Min-Max dapat dilakukan dengan formula berikut:

\[ a'_i = \frac{a_i - min(a)}{max(a)-min(a)}\] Selanjutnya, karena pada R tidak terdapat build-in function dari formula Min-Max ini, maka akan dibuat fungsi manual dengan minMax. Lalu melakukan pengaplikasian fungsi yang sudah dibuat pada data rentan_FCM kita.

# custom function to implement min max scaling
minMax <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}
 
#Normalisasi data menggunakan fungsi custom
rentan_scale <- as.data.frame(lapply(rentan_FCM, minMax))
head(rentan_scale)

#>    perempuan KK_perempuan     lansia pengangguran wajib_ktp     banjir
#> 1 0.05295316    0.2576572 0.04693080   0.21383515 0.2149322 0.52542373
#> 2 0.52031300    0.6466472 0.67779574   0.21787096 0.6517199 0.00000000
#> 3 0.12482581    0.2962708 0.09707556   0.20349840 0.2720766 0.93220339
#> 4 0.70446993    0.4727347 0.46795975   0.04551688 0.5282226 0.03389831
#> 5 0.18362097    0.4507673 0.01395345   0.29043789 0.1014482 0.64406780
#> 6 0.27318041    0.1840792 0.15352836   0.42899989 0.3010948 0.25423729
#>     informal        RS  puskesmas positiv_cov     rawat     isoman
#> 1 0.01814115 0.1428571 0.54545455   0.8162135 0.5774648 0.58285714
#> 2 0.39627308 0.2857143 0.09090909   0.3376162 0.3521127 0.09714286
#> 3 0.03286447 0.1428571 0.63636364   0.9082132 0.7676056 0.83428571
#> 4 0.24102604 0.4285714 0.27272727   0.5649890 0.4366197 0.65142857
#> 5 0.26129471 0.1428571 0.63636364   0.5435508 0.5140845 0.36000000
#> 6 0.41400807 0.0000000 0.72727273   0.7237169 0.4436620 0.70857143

Pemilihan nilai k optimum

Pada project ini, akan digunakan metode elbow untuk melakukan pemilihan nilai k (jumlah cluster) optimal. Di mana akan dipilih nilai k pada saat nilai k ditambah, penurunan Total WSS (Whitin Sum of Squares) tidak terlalu drastis (dapat dikatakan sudah melandai). Secara statistik, dalam melakukan visualisasi elbow digunakan function fviz_nbclust() dari package factoextra.

set.seed(123)
 
elbow <- 
fviz_nbclust(
  x = rentan_scale,
  FUNcluster = kmeans,
  method = "wss",
  k.max = 10
) + labs(subtitle = "Metode Elbow")
elbow

Nilai k optimum dari data rentan_FCM didapatkan pada saat k=3

FCM menggunakan fcm()

Membuat object fcm.rentan dengan menggunakan fungsi fcm() dengan jumlah \(k=3\), \(m=2\) atau nilai membership degree.

library(ppclust)
set.seed(123)
fcm.rentan <- fcm(rentan_scale, 3, m=2, dmetric="sqeuclidean", pw = 2, 
                          alginitv="kmpp", alginitu="imembrand", 
                          nstart=1, iter.max=1000, con.val=1e-09)
summary(fcm.rentan)

#> Summary for 'fcm.rentan'
#> 
#> Number of data objects:  44 
#> 
#> Number of clusters:  3 
#> 
#> Crisp clustering vector:
#>  [1] 1 3 1 2 1 1 1 1 3 3 1 1 3 1 3 2 1 3 1 2 2 2 1 1 1 2 3 3 2 2 2 1 2 1 2 1 3 3
#> [39] 3 3 3 3 1 3
#> 
#> Initial cluster prototypes:
#>           perempuan KK_perempuan    lansia pengangguran wajib_ktp     banjir
#> Cluster 1 0.1605746    0.2174588 0.0000000    0.1790955 0.2569760 0.00000000
#> Cluster 2 0.2864723    0.6662773 0.4948316    0.2670284 0.4809782 0.44067797
#> Cluster 3 0.4331654    0.7510184 0.4848500    0.2040897 0.5733735 0.01694915
#>             informal        RS puskesmas positiv_cov     rawat    isoman
#> Cluster 1 0.01501601 0.1428571 0.9090909   0.5534890 0.5140845 0.2857143
#> Cluster 2 0.27119979 0.1428571 0.4545455   0.6869454 0.6056338 0.3600000
#> Cluster 3 1.00000000 0.5714286 0.1818182   0.3716902 0.4366197 0.3371429
#> 
#> Final cluster prototypes:
#>           perempuan KK_perempuan    lansia pengangguran wajib_ktp     banjir
#> Cluster 1 0.3152556    0.4134893 0.3212723    0.2128704 0.3783466 0.28769114
#> Cluster 2 0.3221162    0.4199689 0.3367700    0.2082273 0.3898127 0.25317872
#> Cluster 3 0.3505734    0.7027618 0.6751806    0.2755074 0.5942018 0.07892542
#>            informal        RS puskesmas positiv_cov     rawat    isoman
#> Cluster 1 0.2153258 0.2756352 0.5032064   0.6137991 0.5298518 0.4764718
#> Cluster 2 0.2287537 0.2644816 0.4857593   0.5844553 0.4991475 0.4445524
#> Cluster 3 0.5672056 0.2309419 0.2612857   0.3440010 0.3417834 0.2199829
#> 
#> Distance between the final cluster prototypes
#>             Cluster 1   Cluster 2
#> Cluster 2 0.005105141            
#> Cluster 3 0.662566819 0.571045341
#> 
#> Difference between the initial and final cluster prototypes
#>             perempuan KK_perempuan     lansia pengangguran   wajib_ktp
#> Cluster 1  0.15468107    0.1960305  0.3212723   0.03377495  0.12137066
#> Cluster 2  0.03564389   -0.2463084 -0.1580616  -0.05880108 -0.09116555
#> Cluster 3 -0.08259198   -0.0482566  0.1903306   0.07141765  0.02082833
#>                banjir    informal         RS   puskesmas positiv_cov
#> Cluster 1  0.28769114  0.20030979  0.1327781 -0.40588453  0.06031003
#> Cluster 2 -0.18749924 -0.04244609  0.1216244  0.03121389 -0.10249015
#> Cluster 3  0.06197627 -0.43279442 -0.3404867  0.07946748 -0.02768923
#>                 rawat      isoman
#> Cluster 1  0.01576731  0.19075756
#> Cluster 2 -0.10648633  0.08455244
#> Cluster 3 -0.09483634 -0.11715994
#> 
#> Root Mean Squared Deviations (RMSD): 0.6042024 
#> Mean Absolute Deviation (MAD): 19.81903 
#> 
#> Membership degrees matrix (top and bottom 5 rows): 
#>   Cluster 1 Cluster 2  Cluster 3
#> 1 0.4882022 0.4145029 0.09729491
#> 2 0.1106537 0.1271162 0.76223009
#> 3 0.4554046 0.3982856 0.14630985
#> 4 0.3712570 0.4013925 0.22735047
#> 5 0.4585425 0.4192312 0.12222637
#> ...
#>    Cluster 1 Cluster 2 Cluster 3
#> 40 0.1454923 0.1573695 0.6971382
#> 41 0.2380096 0.2644525 0.4975380
#> 42 0.1653044 0.1892569 0.6454387
#> 43 0.4238533 0.3844937 0.1916530
#> 44 0.2839744 0.3216587 0.3943669
#> 
#> Descriptive statistics for the membership degrees by clusters
#>           Size       Min        Q1      Mean    Median        Q3       Max
#> Cluster 1   18 0.3910770 0.4354837 0.4579773 0.4508280 0.4640693 0.5855688
#> Cluster 2   11 0.3483277 0.3929432 0.4256993 0.4370806 0.4701132 0.4801915
#> Cluster 3   15 0.3621822 0.4428936 0.5755343 0.5967759 0.7010497 0.7622301
#> 
#> Dunn's Fuzziness Coefficients:
#> dunn_coeff normalized 
#>  0.4115541  0.1173312 
#> 
#> Within cluster sum of squares by cluster:
#>        1        2        3 
#> 6.638132 4.676674 5.472873 
#> (between_SS / total_SS =  27.07%) 
#> 
#> Available components: 
#>  [1] "u"          "v"          "v0"         "d"          "x"         
#>  [6] "cluster"    "csize"      "sumsqrs"    "k"          "m"         
#> [11] "iter"       "best.start" "func.val"   "comp.time"  "inpargs"   
#> [16] "algorithm"  "call"

Untuk mengetahui jumlah kecamatan untuk setiap clusternya digunakan inspeksi pada object fcm$csize

fcm.rentan$csize

#>  1  2  3 
#> 18 11 15

Berikut merupakan nilai untuk masing-masing kecamatan dan juga hasil cluster setiap kecamatannya. Hasil cluster dapat diakses dengan melakukan inspeksi pada object fcm$cluster.

fcm_keep<-data.frame(rentan_scale, fcm.rentan$cluster)
rmarkdown::paged_table(fcm_keep)

• Kluster 1 berisi 18 Kecamatan, mencakup Kecamatan Cakung, Cengkareng, Cilincing, Cipayung, Ciracas, Duren Sawit, Jagakarsa, Jatinegara, Kalideres, Kebonjeruk, Kemayoran, Koja, Kramatjati, Makasar, Pasar Minggu, Penjaringan, Pulogadung, Tanjung Priok.

• Kluster 2 berisi 11 Kecamatan, mencakup Kecamatan Cilandak, Kebayoran Lama, Kembangan, Kepulauan Seribu Selatan, Kepulauan Seribu Utara, Mampang Prapatan, Pademangan, Palmerah, Pancoran, Pasarrebo, dan Pesanggrahan.

• Kluster 3 berisi 15 Kecamatan, mencakup Kecamatan Cempaka Putih, Gambir, Grogolpetamburan, Johar Baru, Kebayoran Baru, Kelapa Gading, Matraman, Menteng, Sawah Besar, Senen, Setia Budi, Tamansari, Tambora, Tanahabang, dan Tebet.

###Interpretasi: Cluster Profiling

Untuk dapat melakukan cluster profiling, dibutuhkan informasi pusat cluster (centroid) masing-masing variabel terhadap cluster yang dihasilkan. Dalam object fcm.rentan, kita dapat melakukan inspeksi menggunakan fcm.rentan$v yang mana merupakan matriks numerin yang berisi pusat masing-masing cluster (centroid).

v <- data.frame(fcm.rentan$v)
rmarkdown::paged_table(v)

Profiling:

Pada klaster 1 ini dapat dilihat bahwa kerentanan tertinggi didapat oleh RS, puskesmas, positiv_cov, rawat, isoman yang ditunjukkan jumlah pasien positif dan dirawat akibat COVID-19 menjadi masalah pada klaster ini, dimana dari data dapat dilihat kecamatan tersebut memiliki jumlah positif dan dirawat akibat covid cukup tinggi. Selanjutnya, banjir juga merupakan masalah yang ada pada cluster ini, di mana beberapa kecamatan yang berada pada cluster 1 ini merupakan daerah yang rentan banjir.

Pada klaster 2 dapat sama-sama diketahui bahwa di antara 3 klaster lainnya, centroid untuk klaster 2 ini cenderung sedang untuk setiap variabelnya. Tidak ada nilai maksimal untuk variabel-variabel pembentuk kerentanan sosial.

Sementara itu pada klaster 3 dapat dilihat bahwa kerentanan tertinggi terdapat pada variabel-variabel perempuan, KK_perempuan, lansia, pengangguran, jumlah penduduk wajib_ktp, serta penduduk yang bekerja di sektor informal. Dapat dilihat dari data pada kecamatan ini rata-rata persentase lansia dan kepala keluarga perempuan yang lebih tinggi dibandingkan klaster lain, yaitu > 9% untuk persentase lansia dan > 21% untuk persentase KK_Perempuan. Pada klaster ini, bagi lansia akan lebih sulit pulih diakibatkan terdapat batasan pergerakan dan ketergantungan terhadap orang lain. Sedangkan pada kepala keluarga perempuan terdapat keterbatasan pilihan pekerjaan, upah yang rendah dan tanggung jawab mengurus rumah tangga.

Berikut merupakan hasil visualisasi kluster dalam plot 2 dimensi.

fcm.rentan2 <- ppclust2(fcm.rentan, "kmeans")
factoextra::fviz_cluster(fcm.rentan2, data = rentan_scale, 
  ellipse.type = "convex",
  palette = "jco",
  repel = TRUE)