Data

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.

Packages

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

library(ggplot2)
library(factoextra)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.2.3

Input Data

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

Gabungan data

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

Hierarchical Clustering

Penggunaan Jarak Mahalanobis

mahala <- as.dist(mahalanobis.dist(lingkungan))
fviz_dist(mahala)

Pemilihan Model Hirarki

#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.

Pemilihan Banyak Cluster

#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.

Clustering

#Interpretasi, utk melihat anggota dari setiap gerombol

hc.data <- eclust(lingkungan, FUNcluster = "hclust", k=3, hc_method = "centroid", hc_metric = "euclidean", graph = F)

Cluster vs order

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.

Barplot Frekuensi Cluster

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

Cluster Plot

fviz_cluster(hc.data)

Dari cluster plot di atas, dapat terlihat bahwa:

  1. 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.

  2. Data masing-masing cluster tidak terlalu menumpuk satu sama lain. artinya, metode clustering telah baik untuk mengelompokkan setiap data.

Sum Square Within (Error)

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.

Non-Hierarchical Clustering (K-Means)

Pemilihan banyaknya cluster

#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.

Clustering

#Ditetapkan k=2
kmeans.data <- eclust(lingkungan, FUNcluster = "kmeans", k=2, graph = F)

Cluster vs Order

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.

Bar Frekuensi Cluster

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.

Cluster Plot

fviz_cluster(kmeans.data)

Dari cluster plot di atas, dapat terlihat bahwa:

  1. 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.

  2. Data masing-masing cluster tidak terlalu menumpuk satu sama lain. artinya, metode clustering telah baik untuk mengelompokkan setiap data.

Sum Square Within (Error)

ssw_k <- kmeans.data$tot.withinss
ssw_k
## [1] 5818.586

Metode Clustering Terbaik

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.

Interpretasi Hasil Clustering

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.