Data yang digunakan adalah data provinsi di Indonesia yang terdiri atas empat peubah, yaitu Indeks Kualitas Udara (IKU), Indeks Kualitas Air (IKA), Indeks Kualitas Air Laut (IKAL), dan Indeks Kualitas Lahan (IKL) pada setiap provinsi. Data ini tersedia pada Laporan Kinerja Direktorat Jenderal Pengendalian Pencemaran dan Kerusakan Lingkungan Tahun 2021. Data ini terdiri dari 34 observasi di setiap peubahnya.
lapply(c("readxl","readr","rvest","factoextra","hrbrthemes",
"gridExtra","StatMatch","scales","dplyr","ggcorrplot",
"glmnet","caret"),
library,character.only=T)[[1]]
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'rvest' was built under R version 4.2.3
##
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
##
## guess_encoding
## Warning: package 'factoextra' was built under R version 4.2.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.2.3
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
## Warning: package 'hrbrthemes' was built under R version 4.2.3
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
## Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
## if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
## Warning: package 'StatMatch' was built under R version 4.2.3
## Loading required package: proxy
##
## Attaching package: 'proxy'
## The following objects are masked from 'package:stats':
##
## as.dist, dist
## The following object is masked from 'package:base':
##
## as.matrix
## Loading required package: survey
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
##
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
##
## dotchart
## Loading required package: lpSolve
## Warning: package 'lpSolve' was built under R version 4.2.3
##
## Attaching package: 'scales'
## The following object is masked from 'package:readr':
##
## col_factor
## Warning: package 'dplyr' was built under R version 4.2.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Warning: package 'ggcorrplot' was built under R version 4.2.3
## Warning: package 'glmnet' was built under R version 4.2.3
## Loaded glmnet 4.1-8
## Warning: package 'caret' was built under R version 4.2.3
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:survival':
##
## cluster
## [1] "readxl" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
library(ggplot2)
library(factoextra)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.2.3
library(rio)
datapsd <- import("https://raw.githubusercontent.com/aidara11/pengantar-sains-data/main/TUGAS%20AKHIR/DATA%20PSD%20UAS.csv")
datapsd
## Provinsi IKU IKA IKAL IKL
## 1 Aceh 89.63 57.14 76.51 76.52
## 2 Sumatera Utara 89.55 53.72 81.43 48.89
## 3 Sumatera Barat 90.22 52.55 83.75 66.24
## 4 Riau 90.13 52.25 77.73 50.22
## 5 Jambi 87.08 48.96 83.58 51.47
## 6 Sumatera Selatan 86.28 58.25 75.53 41.25
## 7 Bengkulu 90.81 49.81 83.61 55.52
## 8 Lampung 85.46 57.77 79.56 33.54
## 9 Bangka Belitung 90.39 58.37 82.71 40.10
## 10 Kepulauan Riau 90.91 55.15 75.68 60.39
## 11 DKI Jakarta 66.52 44.19 75.18 26.25
## 12 Jawa Barat 79.34 43.09 87.42 40.78
## 13 Jawa Tengah 84.60 47.94 83.23 41.51
## 14 DI Yogyakarta 88.59 45.73 83.35 29.66
## 15 Jawa Timur 83.20 53.57 82.46 47.36
## 16 Banten 74.14 54.95 85.92 39.21
## 17 Bali 89.28 54.29 85.14 42.11
## 18 Nusa Tenggara Barat 88.52 45.10 80.22 65.59
## 19 Nusa Tenggara Timur 90.51 58.28 87.36 58.65
## 20 Kalimantan Barat 90.71 54.35 77.83 59.35
## 21 Kalimantan Tengah 90.39 55.34 76.52 75.43
## 22 Kalimantan Selatan 89.15 54.75 76.45 50.26
## 23 Kalimantan Timur 88.84 51.92 85.40 82.21
## 24 Kalimantan Utara 93.43 57.34 81.52 99.96
## 25 Sulawesi Utara 91.27 49.69 82.65 61.94
## 26 Sulawesi Tengah 91.33 55.84 87.36 83.10
## 27 Sulawesi Selatan 89.13 56.82 84.82 55.54
## 28 Sulawesi Tenggara 90.89 53.26 81.60 74.34
## 29 Gorontalo 93.96 53.46 84.80 79.21
## 30 Sulawesi Barat 90.97 56.04 81.52 72.66
## 31 Maluku 90.70 55.56 86.07 90.21
## 32 Maluku Utara 91.64 53.08 87.55 86.58
## 33 Papua Barat 95.60 54.44 81.12 100.00
## 34 Papua 94.02 57.83 70.34 100.00
lingkungan <- data.frame(datapsd$IKU, datapsd$IKA, datapsd$IKAL, datapsd$IKL)
rownames(lingkungan) <- datapsd$Provinsi
lingkungan
## datapsd.IKU datapsd.IKA datapsd.IKAL datapsd.IKL
## Aceh 89.63 57.14 76.51 76.52
## Sumatera Utara 89.55 53.72 81.43 48.89
## Sumatera Barat 90.22 52.55 83.75 66.24
## Riau 90.13 52.25 77.73 50.22
## Jambi 87.08 48.96 83.58 51.47
## Sumatera Selatan 86.28 58.25 75.53 41.25
## Bengkulu 90.81 49.81 83.61 55.52
## Lampung 85.46 57.77 79.56 33.54
## Bangka Belitung 90.39 58.37 82.71 40.10
## Kepulauan Riau 90.91 55.15 75.68 60.39
## DKI Jakarta 66.52 44.19 75.18 26.25
## Jawa Barat 79.34 43.09 87.42 40.78
## Jawa Tengah 84.60 47.94 83.23 41.51
## DI Yogyakarta 88.59 45.73 83.35 29.66
## Jawa Timur 83.20 53.57 82.46 47.36
## Banten 74.14 54.95 85.92 39.21
## Bali 89.28 54.29 85.14 42.11
## Nusa Tenggara Barat 88.52 45.10 80.22 65.59
## Nusa Tenggara Timur 90.51 58.28 87.36 58.65
## Kalimantan Barat 90.71 54.35 77.83 59.35
## Kalimantan Tengah 90.39 55.34 76.52 75.43
## Kalimantan Selatan 89.15 54.75 76.45 50.26
## Kalimantan Timur 88.84 51.92 85.40 82.21
## Kalimantan Utara 93.43 57.34 81.52 99.96
## Sulawesi Utara 91.27 49.69 82.65 61.94
## Sulawesi Tengah 91.33 55.84 87.36 83.10
## Sulawesi Selatan 89.13 56.82 84.82 55.54
## Sulawesi Tenggara 90.89 53.26 81.60 74.34
## Gorontalo 93.96 53.46 84.80 79.21
## Sulawesi Barat 90.97 56.04 81.52 72.66
## Maluku 90.70 55.56 86.07 90.21
## Maluku Utara 91.64 53.08 87.55 86.58
## Papua Barat 95.60 54.44 81.12 100.00
## Papua 94.02 57.83 70.34 100.00
mahala <- as.dist(mahalanobis.dist(lingkungan))
fviz_dist(mahala)
#Menampilkan Dendogram
cmp <- fviz_dend(hclust(dist(lingkungan, method = "euclidean"), method = "complete")) + ggtitle("Complete Linkage")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
avg <- fviz_dend(hclust(dist(lingkungan, method = "euclidean"), method = "average"))+ ggtitle("Average Linkage")
cen <- fviz_dend(hclust(dist(lingkungan, method = "euclidean"), method = "centroid")) + ggtitle("Centroid Linkage")
sin <- fviz_dend(hclust(dist(lingkungan, method = "euclidean"), method = "single")) + ggtitle("Single Linkage")
ggarrange(cmp, avg, cen, sin,
ncol = 2, nrow = 2)
Dari output di atas, metode yang menghasilkan clustering yang paling terdistribusi secara sama (evenly distributed) adalah metode centroid linkage. Untuk itu, akan dilakukan pemilihan banyak cluster berdasarkan metode clustering centroid linkage.
#Centroid Linkage
fviz_nbclust(lingkungan, FUNcluster = hcut, method = "silhouette", hc_method = "centroid", hc_metric="euclidean")
Berdasarkan plot di atas, terlihat bahwa banyaknya cluster yang optimal untuk hiclustering metode centroid linkage sebesar tiga cluster.
#Interpretasi, utk melihat anggota dari setiap gerombol
hc.data <- eclust(lingkungan, FUNcluster = "hclust", k=3, hc_method = "centroid", hc_metric = "euclidean", graph = F)
plot(1:34, hc.data$cluster, xlab="Order", ylab="Cluster", main="Cluster vs Order")
Output di atas memperlihatkan plot cluster terhadap urutan data. Dapat terlihat bahwa data dengan urutan 11 cenderung masuk ke cluster 3, sedangkan data urutan awal masuk ke cluster 2, dan data urutan akhir (21-34) cenderung masuk ke Cluster 1.
library(wesanderson)
## Warning: package 'wesanderson' was built under R version 4.2.3
tab_h <- as.data.frame(table(hc.data$cluster))
tab_h <- tab_h[order(tab_h$Freq),]
barplot(tab_h$Freq, main="Barplot Frekuensi Cluster", ylab="Frekuensi", xlab="Cluster", col=wes_palette(n=3, name="BottleRocket1"), names.arg = tab_h$Var1)
Dari barplot di atas, dapat terlihat bahwa banyaknya frekuensi untuk setiap cluster tidak terlalu jomplang. Artinya, segmentasi data sudah cukup baik
fviz_cluster(hc.data)
Dari cluster plot di atas, dapat terlihat bahwa:
Dim1 dapat menjelaskan 50.1% keragaman dan Dim2 dapat menjelaskan 26.2% keragaman. Artinya, informasi yang diberikan oleh cluster plot di atas sebesar 76.3% dari keseluruhan informasi yang terkandung dalam data.
Data masing-masing cluster tidak terlalu menumpuk satu sama lain. artinya, metode clustering telah baik untuk mengelompokkan setiap data.
Nilai sum square within (Error) akan berguna untuk mendapatkan perbandingan metode clustering terbaik.
X.ss1 <- aggregate(lingkungan, by=list(hc.data$cluster), function(x) sum(scale(x, scale=FALSE)^2))
SS_cluster <- rowSums(X.ss1[-1])
ssw_h <- sum(SS_cluster)
ssw_h
## [1] 4779.705
Nilai di atas merupakan total nilai SSW/Error untuk masing-masing cluster pada penggerombolan berhirarki.
#Penentuan k dengan within sum square
fviz_nbclust(lingkungan, FUNcluster = kmeans, method = "silhouette")
Berdasarkan plot di atas, terlihat bahwa banyaknya cluster yang optimal untuk hclustering metode kmeans sebesar dua cluster.
#Ditetapkan k=2
kmeans.data <- eclust(lingkungan, FUNcluster = "kmeans", k=2, graph = F)
plot(1:34, kmeans.data$cluster, xlab="Order", ylab="Cluster", main="Cluster vs Order")
Output di atas memperlihatkan plot cluster terhadap urutan data. Dapat terlihat bahwa data dengan urutan awal cenderung masuk ke cluster ke-2, sedangkan data urutan akhir masuk ke cluster ke-1.
tab_k <- as.data.frame(table(kmeans.data$cluster))
tab_kO <- tab_k[order(tab_k$Freq),]
barplot(tab_kO$Freq, main="Barplot Frekuensi Cluster", ylab="Frekuensi", xlab="Cluster", col=wes_palette(n=2, name="BottleRocket1"), names.arg = tab_kO$Var1)
Dari barplot di atas, dapat terlihat bahwa banyaknya frekuensi untuk setiap cluster tidak terlalu jomplang. Artinya, segmentasi data sudah cukup baik.
fviz_cluster(kmeans.data)
Dari cluster plot di atas, dapat terlihat bahwa:
Dim1 dapat menjelaskan 50.1% keragaman dan Dim2 dapat menjelaskan 26.2% keragaman. Artinya, informasi yang diberikan oleh cluster plot di atas sebesar 76.3% dari keseluruhan informasi yang terkandung dalam data.
Data masing-masing cluster tidak terlalu menumpuk satu sama lain. artinya, metode clustering telah baik untuk mengelompokkan setiap data.
ssw_k <- kmeans.data$tot.withinss
ssw_k
## [1] 5818.586
sst <- sum(lingkungan^2)
ssb_h <- sst-ssw_h
ssb_k <- sst-ssw_k
Fh <- (ssb_h/3)/(ssw_h/30) #34-1-3
Fk <- (ssb_k/2)/(ssw_k/31) #34-1-2
best_method <- data.frame("Hierarchical"=Fh, "KMeans"=Fk)
rownames(best_method)="Nilai Kebaikan Model"
best_method
## Hierarchical KMeans
## Nilai Kebaikan Model 1524.684 1938.544
Dari nilai di atas, dapat terlihat bahwa nilai kebaikan model K-means setelah distandarisasi lebih besar daripada nilai kebaikan model Hirarki. Artinya, metode Kmeans dipilih untuk penggerombolan data.
df_kmeans <- aggregate(lingkungan, by=list(cluster=kmeans.data$cluster), FUN = mean)
cbind(df_kmeans,"Frekuensi"=tab_k$Freq)
## cluster datapsd.IKU datapsd.IKA datapsd.IKAL datapsd.IKL Frekuensi
## 1 1 91.66308 54.90769 81.85077 83.57385 13
## 2 2 86.45571 52.23952 81.51714 47.59952 21
Cluster 1 merupakan segmentasi provinsi dengan indeks kualitas udara tinggi, indeks kualitas air di atas rata-rata, indeks kualitas air laut rata-rata, dan indeks kualitas lahan tinggi. Artinya, segmentasi ini berisi provinsi-provinsi dengan indeks kualitas lingkungan hidup yang sudah baik.
Cluster 2 merupakan segmentasi provinsi dengan indeks kualitas udara rendah, indeks kualitas air di bawah rata-rata, indeks kualitas air laut rata-rata, dan indeks kualitas lahan rendah. Artinya, segmentasi ini berisi provinsi-provinsi dengan indeks kualitas lingkungan hidup yang buruk sehingga perlu penanganan lebih lanjut.