Library:

> # install.packages("knitr")
> # install.packages("rmarkdown")
> # install.packages("prettydoc")
> # install.packages("equatiomatic")
> library(knitr)
> library(rmarkdown)
> library(prettydoc)
> library(equatiomatic)

1 PENDAHULUAN

1.1 Latar Belakang

Indeks Pembangunan Manusia (IPM) merupakan indikator komposit tunggal yang mengukur tiga dimensi pokok pembangunan manusia yang dinilai mampu mencerminkan kemampuan dasar (basic capabilities) penduduk (Setiawan & Hakim, 2013). Ketiga kemampuan dasar itu adalah umur panjang dan sehat, berpengetahuan dan berketerampilan, serta akses terhadap sumber daya yang dibutuhkan untuk mencapai standar hidup layak. Ketiga dimensi tersebut memiliki pengertian sangat luas karena terkait banyak faktor. Setiap indikator komponen penghitungan IPM dapat dimanfaatkan untuk mengukur capaian pembangunan kualitas hidup manusia (BPS, 2021).

Dimensi umur panjang dan hidup sehat diukur dengan indikator umur harapan hidup saat lahir. Dimensi pengetahuan diukur dengan gabungan indikator harapan lama sekolah dan rata-rata lama sekolah. Adapun untuk mengukur dimensi standar hidup layak digunakan indikator kemampuan daya beli. Kemampuan daya beli masyarakat terhadap sejumlah kebutuhan pokok yang dilihat dari rata-rata besarnya pengeluaran rill per kapita sebagai pendekatan pendapatan yang mewakili capaian pembangunan untuk standar hidup layak. Pengelompokan wilayah berdasarkan indikator-indikator ini dapat membantu dalam mengidentifikasi wilayah yang belum mencapai tingkat pembangunan optimal.

Semakin tinggi angka IPM suatu wilayah, semakin mecerminkan keberhasilan wilayah tersebut dalam proses pembangunan. Untuk meningkatkan keberhasilan ini, diperlukan identifikasi terhadap wilayah-wilayah yang membutuhkan perbaikan. Salah satu langkah yang dapat dilakukan adalah mengelompokkan wilayah-wilayah yang belum mencapai tingkat pembangunan yang optimal. Pengelompokan ini dapat membantu pemerintah untuk menyesuaikan program pembangunan sesuai dengan karakteristik dan kebutuhan masing-masing wilayah, yang pada akhirnya akan meningkatkan efektivitas kebijakan dan pemerataan pembangunan di seluruh wilayah.

1.2 Latar Belakang Metode

Indeks Pembangunan Manusia (IPM) merupakan indikator penting dalam pengukuran tingkat keberhasilan pembangunan kualitas hidup manusia (Sikana & Wijayanto, 2021). Pengelompokan Indeks Pembangunan Manusia (IPM) bertujuan untuk membagi wilayahwilayah ke dalam beberapa kelompok yang memiliki tingkat keseragaman karakteristik IPM yang tinggi dalam setiap kelompok dan keberagaman karakteristik IPM yang tinggi antar kelompok. Karakteristik wilayah yang kompleks dan beragam sering kali menjadi tantangan, sehingga dibutuhkan metode analisis yang dapat mengelompokkan wilayah secara objektif, terukur, dan efisien.

Salah satu metode yang cocok untuk tujuan tersebut adalah analisis cluster K-Means. Metode ini mampu mengelompokkan wilayah berdasarkan kemiripan nilai indikator IPM, seperti umur harapan hidup, harapan lama sekolah, rata-rata lama sekolah, dan daya beli. K-Means bekerja dengan mengidentifikasi pola data dan membagi wilayah ke dalam kelompok (cluster) yang memiliki karakteristik serupa.

Dengan analisis cluster K-Means, pemerintah dapat mengidentifikasi wilayah-wilayah yang belum mencapai tingkat pembangunan optimal secara lebih terstruktur. Hasil pengelompokan ini memberikan landasan yang jelas untuk merancang program pembangunan yang lebih tepat sasaran, efektif, dan berorientasi pada pemerataan pembangunan manusia.

1.3 Cluster K-Means

Teknik clustering adalah suatu metode dalam analisa data untuk menyelesaikan persoalan tentang pengelompokkan data (Maori & Evanita, 2023). Salah satu metode yang digunakan adalah K-Means, yang mempunyai kemampuan dalam klasifikasi data yang cukup besar jumlahnya cukup relatif cepat dan efisien. K-means merupakan salah satu metode pengelompokkan data non-hierarki yang mempartisi data yang ada ke dalam bentuk dua atau lebih kelompok (Asroni dkk., 2018).

Terdapat beberapa langkah utama dalam melakukan K-Means clustering, yaitu:

Menentukan jumlah cluster dan menentukan koordinat titik tengah cluster, dalam analisis ini menggunakan metode elbow.
Penentuan nilai dari cluster untuk dijadikan acuan dalam melakukan perhitungan jarak objek ke centroid, perhitungan jarak mengacu pada rumus euclidean.
Dilakukan pengelompokkan centroid sesuai dengan hasil dari jarak antar centroid tersebut.

1.4 Rumusan Masalah

Rumusan masalah dalam penelitian ini adalah sebagai berikut.

Bagaimana kelompok atau cluster terbentuk dari 38 Kabupaten/Kota di Jawa Timur pada tahun 2021 berdasarkan indikator IPM?
Kelompok atau cluster mana yang memiliki tingkat IPM lebih tinggi, sedang, dan rendah?

1.5 Tujuan Penelitian

Berdasarkan masalah yang telah dirumuskan, tujuan penelitian ini antara lain:

Menganalisis dan menentukan kelompok atau cluster yang terbentuk dari 38 Kabupaten/Kota di Jawa Timur pada tahun 2021 berdasarkan indikator Indeks Pembangunan Manusia (IPM) menggunakan metode K-Means Clustering.
Menentukan kelompok wilayah yang memiliki tingkat pembangunan manusia lebih tinggi dan lebih rendah.

1.6 Batasan Masalah

Batasan masalah yang digunakan pada penelitian ini adalah sebagai berikut.

Data yang digunakan adalah data sekunder Indeks Pembangunan Manusia pada kabupaten/kota di Jawa Timur tahun 2021 yang didapatkan dari situs resmi Badan Pusat Statistik (BPS) Jawa Timur.

1.7 Data

Data yang digunakan berasal dari data 38 Kabupaten/Kota di Jawa Timur pada tahun 2021 yang diperoleh dari situs resmi Badan Pusat Statistik (BPS) Jawa Timur yang meliputi sejumlah indikator Indeks Pembangunan Manusia (IPM), yakni Angka Harapan Hidup (AHH), Harapan Lama Sekolah (HLS), Pengeluaran per Kapita (PK), Rata-rata Lama Sekolah (RLS) dan Pembangunan.

Adapun data tersebut dapat dilihat pada tabel berikut:

Indikator Indeks Pembangunan Manusia (IPM) terdiri dari 5 indikator, yaitu:

\(X_1\) : Angka Harapan Hidup (Tahun)
\(X_2\) : Harapan Lama Sekolah (Tahun)
\(X_3\) : Pengeluaran per Kapita (Rupiah)
\(X_4\) : Rata-rata Lama Sekolah (Tahun)
\(X_5\) : Pembangunan (%)

2 SOURCE CODE

2.1 Library

Digunakan untuk membaca file excel dengan format file spreadsheet Microsoft Excel (xlsx).

> library(readxl)

Digunakan untuk mengevaluasi kualitas clustering, dengan menghitung berbagai metrik validasi internal seperti konektivitas, Dunn index, dan Silhouette index, untuk membantu menentukan jumlah cluster yang optimal dan mengevaluasi seberapa baik data dikelompokkan.

> library(clValid)

Digunakan dalam analisis K-Means clustering untuk membantu dalam visualisasi hasil klastering dan menentukan jumlah klaster optimal.

> library(factoextra)

Digunakan untuk membuat tabel interaktif yang mudah digunakan dalam analisis data

> library(DT)

2.2 Import Data

Mengimpor data dari file excel dengan cara membaca file Excel menggunakan fungsi read_excel(), lalu menyimpannya dalam “data_cluster”

> data_cluster <- read_excel("D:/LAPRAK NAZWA/Data IPM Jawa Timur 2021.xlsx")
> data_cluster
# A tibble: 38 × 6
   `Kabupaten/Kota`        AHH   HLS    PK   RLS Pertumbuhan
   <chr>                 <dbl> <dbl> <dbl> <dbl>       <dbl>
 1 Kabupaten Pacitan      72.1  12.6  8887  7.61        0.26
 2 Kabupaten Ponorogo     72.8  13.7  9851  7.55        0.35
 3 Kabupaten Trenggalek   73.9  12.5  9743  7.56        0.46
 4 Kabupaten Tulungagung  74.2  13.3 10807  8.34        0.21
 5 Kabupaten Blitar       73.6  12.6 10757  7.5         0.67
 6 Kabupaten Kediri       72.6  13.4 11127  8.08        0.71
 7 Kabupaten Malang       72.6  13.2 10163  7.43        0.34
 8 Kabupaten Lumajang     70.2  11.9  9203  6.67        0.93
 9 Kabupaten Jember       69.3  13.4  9410  6.49        0.31
10 Kabupaten Banyuwangi   70.7  13.1 12217  7.42        1.08
# ℹ 28 more rows
> View(data_cluster)

2.3 Stuktur Data

Digunakan untuk menampilkan struktur objek data seperti tipe data objek, jumlah elemen dalam objek, dan informasi lainnya. Struktur data digunakan untuk memastikan tipe data yang akan dianalisis sudah benar.

str(data_cluster)

2.4 Memeriksa Missing Value (NA)

Digunakan untuk memeriksa apakah terdapat nilai NA (missing values) pada kolom ke-2 hingga ke-6 dalam dataset, yang akan mengembalikan TRUE jika ada NA, dan FALSE jika tidak ada NA.

data_cluster[,2:6] %>%
anyNA()

2.5 Memeriksa Jumlah NA pada Setiap Kolom

Digunakan untuk menghitung jumlah nilai NA (missing values) pada setiap kolom dari kolom ke-2 hingga ke-6 dalam dataset.

data_cluster[,2:6] %>%
is.na() %>%
colSums()

2.6 Statistika Deskriptif

Digunakan untuk memberikan pemahaman awal mengenai data yang akan diuji. Menampilkan statistik deskriptif (informasi seperti mean, median, maksimum, minimum, dan kuartil) dari tiap variabel yang ada.

summary(data_cluster)

2.7 Menentukan Cluster Optimal (Metode Elbow)

Digunakan untuk menentukan jumlah cluster optimal dalam K-Means clustering dengan metode Elbow. Fungsi “fviz_nbclust()” dari library factoextra dipakai untuk memvisualisasikan jumlah cluster dalam dataset.

fviz_nbclust(data_cluster[,2:6], kmeans, method = "wss") + labs(title = "Elbow Method")

2.8 Clustering K Means

Digunakan fungsi eclust() untuk menjalankan clustering K-Means dari dataset dengan jumlah cluster sebanyak tiga. Setelah hasil clustering diperoleh, pipe operator %>% digunakan untuk mengakses objek hasil clustering dan mengubah nama kolom pada matriks centers yang berisi nilai rata-rata dari masing-masing cluster.

K_Means = eclust(data_cluster[, 2:6], "kmeans", 3, graph = F) %>%
  {colnames(.$centers) <- c("AHH", "HLS", "PK", "RLS", "Pertumbuhan"); .}
  
K_Means

2.9 Cluster Plot

fviz_cluster(K_Means, 
             data = data_cluster[, 2:6], 
             palette = c("#E41A1C", "#377EB8", "#4DAF4A"),
             geom = c("point", "text"), 
             main = "Plot Clustering K-Means")

2.10 Indeks Validitas

Setelah melakukan clustering menggunakan K-Means, penting untuk mengevaluasi kualitas hasil clustering dengan menggunakan indeks validitas, karena membantu menilai seberapa baik data telah dikelompokkan dengan mengukur seberapa erat dan terpisahnya setiap cluster yang terbentuk. Indeks ini memastikan bahwa data dalam satu cluster memiliki kesamaan yang tinggi, sementara data antar cluster memiliki perbedaan yang signifikan.

cl_valid_result <- clValid(data.frame(data_cluster[, 2:6]), nClust = 3, clMethods = "kmeans", 
                           validation = "internal", metric = "euclidean", method = "complete")
summary(cl_valid_result)

3 HASIL DAN INTERPRETASI

3.1 Stuktur Data

> str(data_cluster)
tibble [38 × 6] (S3: tbl_df/tbl/data.frame)
 $ Kabupaten/Kota: chr [1:38] "Kabupaten Pacitan" "Kabupaten Ponorogo" "Kabupaten Trenggalek" "Kabupaten Tulungagung" ...
 $ AHH           : num [1:38] 72.1 72.8 73.9 74.2 73.6 ...
 $ HLS           : num [1:38] 12.6 13.7 12.5 13.3 12.6 ...
 $ PK            : num [1:38] 8887 9851 9743 10807 10757 ...
 $ RLS           : num [1:38] 7.61 7.55 7.56 8.34 7.5 8.08 7.43 6.67 6.49 7.42 ...
 $ Pertumbuhan   : num [1:38] 0.26 0.35 0.46 0.21 0.67 0.71 0.34 0.93 0.31 1.08 ...

Interpretasi: Berdasarkan output diatas, dapat dilihat bahwa tipe data pada variabel AHH sampai dengan Pertumbuhan adalah numerik.

3.2 Memeriksa Missing Value (NA)

> data_cluster[,2:6] %>%
+ anyNA()
[1] FALSE

Interpretasi: Setelah dilakukan pemeriksaan terhadap missing value (NA) dalam dataset didapatkan bahwa tidak terdapat missing value.

3.3 Memeriksa Jumlah NA pada Setiap Kolom

> data_cluster[,2:6] %>%
+ is.na() %>%
+ colSums()
        AHH         HLS          PK         RLS Pertumbuhan 
          0           0           0           0           0

Interpretasi: Didapatkan hasil bahwa tidak terdapat missing value (NA) pada seluruh indikator sehingga dapat dianalisis lebih lanjut.

3.4 Statistika Deskriptif

> summary(data_cluster[,2:6])
      AHH             HLS              PK             RLS        
 Min.   :66.89   Min.   :11.73   Min.   : 8673   Min.   : 4.860  
 1st Qu.:70.44   1st Qu.:12.72   1st Qu.:10038   1st Qu.: 7.200  
 Median :72.45   Median :13.32   Median :11260   Median : 7.695  
 Mean   :71.72   Mean   :13.40   Mean   :11569   Mean   : 8.061  
 3rd Qu.:72.84   3rd Qu.:13.77   3rd Qu.:12723   3rd Qu.: 9.220  
 Max.   :74.18   Max.   :15.75   Max.   :17862   Max.   :11.370  
  Pertumbuhan    
 Min.   :0.1000  
 1st Qu.:0.3175  
 Median :0.4750  
 Mean   :0.4937  
 3rd Qu.:0.6675  
 Max.   :1.0800

Interpretasi: Berdasarkan statistika deskriptif di atas, didapatkan informasi:

Kabupaten/Kota yang memiliki angka harapan hidup tertinggi adalah Kota Surabaya sebesar 74.18 tahun dan terendah adalah Kabupaten Bondowoso sebesar 66.89 tahun. Sebagian besar Kabupaten/Kota di Jawa Timur tahun 2021 memiliki angka harapan hidup sebesar 71.72 tahun.
Kabupaten/Kota yang memiliki harapan lama sekolah tertinggi adalah Kota Malang sebesar 15.75 tahun dan terendah adalah Kabupaten Bangkalan sebesar 11.73 tahun. Sebagian besar Kabupaten/Kota di Jawa Timur tahun 2021 memiliki harapan lama sekolah sebesar 13.40 tahun.
Kabupaten/Kota yang memiliki pengeluaran per kapita tertinggi adalah Kota Surabaya sebesar Rp 17.862 dan terendah adalah Kabupaten Bangkalan sebesar Rp 8.673. Sebagian besar Kabupaten/Kota di Jawa Timur tahun 2021 memiliki pengeluaran per kapita sebesar Rp 11.569.
Kabupaten/Kota yang memiliki rata-rata lama sekolah tertinggi adalah Kota Madiun sebesar 11.37 tahun dan terendah adalah Kabupaten Sampang sebesar 4.86 tahun. Sebagian besar Kabupaten/Kota di Jawa Timur tahun 2021 memiliki rata-rata lama sekolah sebesar 8.06 tahun.
Kabupaten/Kota yang memiliki pertumbuhan tertinggi adalah Kabupaten Banyuwangi sebesar 1.08 persen dan terendah adalah Kota Surabaya sebesar 0.1 persen. Sebagian besar Kabupaten/Kota di Jawa Timur tahun 2021 memiliki pertumbuhan sebesar 0.4937 persen.

3.5 Menentukan Cluster Optimal (Metode Elbow)

> fviz_nbclust(data_cluster[,2:6], kmeans, method = "wss") + labs(title = "Elbow Method")

Interpretasi: Berdasarkan plot diatas, dapat dilihat bahwa titik yang menurun tajam dibandingkan dengan titik yang lain adalah titik ketiga dan setelah titik ini pengurangan TWSS tidak lagi signifikan, sehingga dengan menggunakan metode wss (elbow) dapat diperoleh banyaknya cluster adalah 3 cluster.

3.6 Clustering K Means

> K_Means = eclust(data_cluster[, 2:6], "kmeans", 3, graph = F) %>%
+   {colnames(.$centers) <- c("AHH", "HLS", "PK", "RLS", "Pertumbuhan"); .}
>   
> K_Means
K-means clustering with 3 clusters of sizes 16, 4, 18

Cluster means:
       AHH      HLS        PK       RLS Pertumbuhan
1 72.35375 13.65687 12360.312  8.754375   0.5443750
2 73.60750 14.97750 16299.500 10.750000   0.4225000
3 70.73278 12.81889  9813.389  6.846667   0.4644444

Clustering vector:
 [1] 3 3 3 3 3 1 3 3 3 1 3 3 3 3 2 1 1 1 1 1 1 3 3 1 1 3 3 3 3 1 1 2 1 1 1 2 2 1

Within cluster sum of squares by cluster:
[1] 10856191  5578924 10286450
 (between_SS / total_SS =  85.3 %)

Available components:

 [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
 [6] "betweenss"    "size"         "iter"         "ifault"       "silinfo"     
[11] "nbclust"      "data"

Interpretasi: Berdasarkan analisis cluster menggunakan K-Means diatas, didapatkan cluster 1 berjumlah 16 Kabupaten/Kota, cluster 2 berjumlah 4 Kabupaten/Kota, cluster 3 berjumlah 18 Kabupaten/Kota

Angka Haraoan Hidup (AHH) tertinggi terdapat pada cluster 2, yang dapat mencerminkan wilayah dengan kualitas hidup atau akses kesehatan yang lebih baik. Harapan Lama Sekolah (HLS) tertinggi pada cluster 2, Pengeluaran per Kapita (PK) tertinggi pada cluster 2, Rata-rata Lama Sekolah (RLS) tertinggi pada cluster 2, dimana Cluster 2 menonjol dalam empat indikator utama (AHH, HLS, PK, dan RLS), menunjukkan bahwa wilayah dalam cluster ini cenderung memiliki kualitas hidup, tingkat pendidikan dan kesejahteraan ekonomi yang lebih tinggi dan baik. Cluster 1 unggul dalam indikator Pertumbuhan, menunjukkan wilayah dengan dinamika perkembangan yang lebih pesat. Cluster 3 memiliki karakteristik yang kurang menarik karena tidak unggul pada salah satu dari lima indikator, namun berdasarkan banyaknya Kabupaten/Kota, cluster 3 memiliki jumlah anggota terbanyak.

3.7 Cluster Plot

> fviz_cluster(K_Means, 
+              data = data_cluster[, 2:6], 
+              palette = c("#E41A1C", "#377EB8", "#4DAF4A"),
+              geom = c("point", "text"), 
+              main = "Plot Clustering K-Means")

Interpretasi: Cluster 1 dan 3 beririsan, menunjukkan bahwa ada beberapa objek dengan karakteristik yang mirip di antara cluster tersebut.

3.8 Indeks Validitas

> cl_valid_result <- clValid(data.frame(data_cluster[, 2:6]), nClust = 3, clMethods = "kmeans", 
+                            validation = "internal", metric = "euclidean", method = "complete")
> summary(cl_valid_result)

Clustering Methods:
 kmeans 

Cluster sizes:
 3 

Validation Measures:
                           3
                            
kmeans Connectivity  12.6123
       Dunn           0.0481
       Silhouette     0.5613

Optimal Scores:

             Score   Method Clusters
Connectivity 12.6123 kmeans 3       
Dunn          0.0481 kmeans 3       
Silhouette    0.5613 kmeans 3

Interpretasi: Clustering cukup baik, tetapi belum optimal karena memiliki Silhouette Score 0.5613. Dunn Index rendah dan Connectivity tinggi mengindikasikan bahwa ada masalah dengan pemisahan antar-cluster dan kepadatan cluster.

4 PENUTUP

4.1 Kesimpulan

Kesimpulan yang diperoleh berdasarkan hasil analisis, adalah sebagai berikut.

Tiga cluster terbentuk berdasarkan indikator Indeks Pembangunan Manusia (IPM):

Cluster 1 (16 Kabupaten/Kota): Wilayah ini memiliki performa unggul dalam indikator pertumbuhan ekonomi, mencerminkan dinamika pembangunan yang cepat.
Cluster 2 (4 Kabupaten/Kota): Wilayah ini memiliki tingkat Harapan Hidup (AHH), Harapan Lama Sekolah (HLS), Pengeluaran per Kapita (PK), dan Rata-rata Lama Sekolah (RLS) tertinggi, mencerminkan kualitas hidup, pendidikan, dan kesejahteraan ekonomi yang lebih baik.
Cluster 3 (18 Kabupaten/Kota): Wilayah dengan jumlah anggota terbanyak, tetapi tidak unggul dalam indikator apa pun, mencerminkan wilayah dengan pembangunan manusia yang relatif biasa saja atau tidak menonjol.

Cluster 2 merupakan cluster tinggi, terdiri dari wilayah-wilayah yang menunjukkan indikator IPM terbaik, termasuk AHH, HLS, PK, dan RLS. Wilayah dalam cluster ini dapat dianggap sebagai pusat kesejahteraan dan kualitas hidup yang tinggi. Cluster 3 merupakan cluster rendah, wilayah dengan karakteristik pembangunan yang kurang menonjol (stagnan), mencerminkan tantangan dalam meningkatkan kualitas pendidikan, kesejahteraan ekonomi, dan kesehatan. Cluster 1 merupakan cluster sedang.

4.2 Saran

Pemerintah perlu meningkatkan akses pendidikan dan kesehatan di wilayah Cluster 3 seperti membangun infrastruktur sekolah dan meningkatkan fasilitas kesehatan. Selain itu, pemberdayaan ekonomi lokal dapat dilakukan melalui dukungan kepada UMKM dan investasi infrastruktur yang sesuai dengan potensi daerah. Hal ini bertujuan untuk meningkatkan kualitas hidup masyarakat sekaligus mengurangi kesenjangan antarwilayah.
Berdasarkan analisis yang sudah dilakukan, dapat dilakukan analisis lanjutan dengan memperbanyak indikator-indikator lain terkait Indeks Pembangunan Manusia.

5 DAFTAR PUSTAKA

Asroni, Fitri, H., Prasetyo, E. (2018). Penerapan Metode Clustering dengan Algoritma K-Means pada Pengelompokkan Data Calon Mahasiswa Baru di Universitas Muhammadiyah Yogyakarta (Studi Kasus: Fakultas Kedokteran dan Ilmu Kesehatan, dan Fakultas Ilmu Sosial dan Ilmu Politik). Semesta Teknika. 21(1), hal 60-63.

BPS. (2021). Indeks Pembangunan Manusia 2021. Jakarta: Badan Pusat Statistik.

Maori, N. A. & Evanita. (2023). METODE ELBOW DALAM OPTIMASI JUMLAH CLUSTER PADA K-MEANS CLUSTERING. Jurnal SIMETRIS. 14(2), hal 278-279.

Setiawan, M. B. & Hakim, A. (2013). INDEKS PEMBANGUNAN MANUSIA INDONESIA. Jurnal Economica. 9(1), hal 20.

Sikana, A. M. & Wijayanto, A. W. (2021). Analisis Perbandingan Pengelompokan Indeks Pembangunan Manusia Indonesia Tahun 2019 dengan Metode Partitioning dan Hierarchical Clustering . Jurnal Ilmu Computer. 14(2), hal 67.

Analisis K-Means untuk Pengelompokkan Kabupaten/Kota di Jawa Timur Berdasarkan Indeks Pembangunan Manusia Tahun 2021

Nazwa Anindya Putri

27 November 2024