R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(psych)
library(ggcorrplot)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ lubridate 1.9.4     ✔ stringr   1.5.1
## ✔ purrr     1.0.2     ✔ tibble    3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%()   masks psych::%+%()
## ✖ ggplot2::alpha() masks psych::alpha()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(scatterplot3d)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
# Baca Data
data_IKK <- read.csv2("C:/Users/Me/OneDrive/Dokumen/IKK.csv")
head(data_IKK)
##          Wilayah  IKK   IPM   TPT   PPK   RLS
## 1 Kab Pandeglang 1.21 65.84  9.24  8827  7.13
## 2      Kab Lebak 1.35 64.71  8.55  8854  6.59
## 3  Kab Tangerang 1.06 72.97  7.88 12427  8.92
## 4     Kab Serang 0.56 67.75 10.61 10916  7.78
## 5 Kota Tangerang 0.72 78.90  7.16 14909 10.84
## 6   Kota Cilegon 0.42 73.95  8.10 13185 10.34
str(data_IKK)
## 'data.frame':    119 obs. of  6 variables:
##  $ Wilayah: chr  "Kab Pandeglang" "Kab Lebak" "Kab Tangerang" "Kab Serang" ...
##  $ IKK    : num  1.21 1.35 1.06 0.56 0.72 0.42 0.77 0.32 2.32 2.11 ...
##  $ IPM    : num  65.8 64.7 73 67.8 78.9 ...
##  $ TPT    : num  9.24 8.55 7.88 10.61 7.16 ...
##  $ PPK    : int  8827 8854 12427 10916 14909 13185 13709 15997 10511 16002 ...
##  $ RLS    : num  7.13 6.59 8.92 7.78 10.84 ...
describe(data_IKK)
##          vars   n     mean      sd   median  trimmed     mad     min      max
## Wilayah*    1 119    60.00   34.50    60.00    60.00   44.48    1.00   119.00
## IKK         2 119     1.42    0.68     1.33     1.36    0.59    0.32     3.72
## IPM         3 119    73.71    5.31    72.97    73.42    5.09   63.39    87.69
## TPT         4 119     6.12    2.28     5.92     6.12    2.49    1.36    10.78
## PPK         5 119 12545.93 3752.60 12379.00 12509.85 3770.25 4173.00 24221.00
## RLS         6 119    10.18    2.67     9.62    10.08    3.29    6.12    15.76
##             range skew kurtosis     se
## Wilayah*   118.00 0.00    -1.23   3.16
## IKK          3.40 0.94     0.77   0.06
## IPM         24.30 0.53    -0.53   0.49
## TPT          9.42 0.03    -0.68   0.21
## PPK      20048.00 0.18    -0.08 344.00
## RLS          9.64 0.26    -1.34   0.24
summary(data_IKK)
##    Wilayah               IKK             IPM             TPT        
##  Length:119         Min.   :0.320   Min.   :63.39   Min.   : 1.360  
##  Class :character   1st Qu.:0.930   1st Qu.:69.78   1st Qu.: 4.495  
##  Mode  :character   Median :1.330   Median :72.97   Median : 5.920  
##                     Mean   :1.425   Mean   :73.71   Mean   : 6.124  
##                     3rd Qu.:1.735   3rd Qu.:76.77   3rd Qu.: 7.830  
##                     Max.   :3.720   Max.   :87.69   Max.   :10.780  
##       PPK             RLS        
##  Min.   : 4173   Min.   : 6.120  
##  1st Qu.: 9974   1st Qu.: 7.755  
##  Median :12379   Median : 9.620  
##  Mean   :12546   Mean   :10.184  
##  3rd Qu.:15105   3rd Qu.:12.710  
##  Max.   :24221   Max.   :15.760
dat <- data_IKK
# Exploratory Data Analysis
data_long <- data_IKK %>%
  pivot_longer(
    cols = c(IKK, IPM, TPT, PPK, RLS),
    names_to = "Variable",
    values_to = "Value"
  )

ggplot(data_long, aes(x = Value, fill = Variable)) +
  geom_histogram(bins = 30, color = "black") +
  facet_wrap(~ Variable, scales = "free") +
  labs(title = "Distribution of Variables",
       x = "Value",
       y = "Count") +
  theme_minimal()

# Standardization
dat_mds <- dat %>% select(IKK, IPM, TPT, PPK, RLS)
data_scaled <- scale(dat_mds)
head(data_scaled)
##             IKK        IPM       TPT        PPK         RLS
## [1,] -0.3177316 -1.4805002 1.3674590 -0.9910274 -1.14343261
## [2,] -0.1103909 -1.6931774 1.0646813 -0.9838324 -1.34563686
## [3,] -0.5398823 -0.1385636 0.7706799 -0.0316934 -0.47316297
## [4,] -1.2803848 -1.1210193 1.9686261 -0.4343472 -0.90003861
## [5,] -1.0434240  0.9775211 0.4547380  0.6297141  0.24578546
## [6,] -1.4877255  0.0458821 0.8672177  0.1702997  0.05855931
# Method 1 : Multidimensional Scaling 
matriks_jarak <- as.matrix(dist(data_scaled))
A <- matriks_jarak^2
n <- nrow(matriks_jarak)
I <- diag(n)
J <- matrix(rep(1, n), nrow = n, ncol = n)
V <- I - (1/n) * J
aa <- V %*% A
BB <- aa %*% V
B <- (-1/2) * BB

eigen_hasil <- eigen(B)
nilai_eigen <- eigen_hasil$values
vektor_eigen <- eigen_hasil$vectors

cumulative_variance <- cumsum(nilai_eigen) / sum(nilai_eigen)
cumulative_variance
##   [1] 0.4483767 0.6859893 0.8517372 0.9580564 1.0000000 1.0000000 1.0000000
##   [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [15] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [22] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [29] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [36] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [43] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [50] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [57] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [64] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [71] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [78] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [85] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [92] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##  [99] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [106] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [113] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
fit <- cmdscale(matriks_jarak, k = 3)
colnames(fit) <- c("MDS1", "MDS2", "MDS3")
dat$MDS1 <- fit[, 1]
dat$MDS2 <- fit[, 2]
dat$MDS3 <- fit[, 3]
#STRESS untuk setiap dimensi
max_dim <- 5
stress_values <- numeric(max_dim)
for (k in 1:max_dim) {
  fit_k <- cmdscale(matriks_jarak, k = k)
  disparities_k <- as.matrix(dist(fit_k))
  stress_values[k] <- sqrt(sum((matriks_jarak - disparities_k)^2) / sum(matriks_jarak^2))
}
stress_values
## [1] 4.810098e-01 2.701668e-01 1.332582e-01 4.216117e-02 8.407824e-16
plot(1:max_dim, stress_values, type = "b", pch = 19, col = "blue",
     xlab = "Dimensi", ylab = "Nilai STRESS",
     main = "Grafik STRESS setiap Dimensi",
     ylim = c(min(stress_values) - 0.01, max(stress_values) + 0.01))
abline(h = 0.01, col = "red", lty = 2)

#Visualisasi MDS 2D
ggplot(dat, aes(x = MDS1, y = MDS2, label = Wilayah)) +
  geom_point(size = 3, color = "blue") +
  geom_text(vjust = -1, size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Visualisasi MDS 2D Kabupaten/Kota Provinsi Jawa Tengah berdasarkan IKK 2023",
    x = "Dimensi 1",
    y = "Dimensi 2"
  ) +
  theme_minimal()

stress_2d <- sqrt(sum((matriks_jarak - as.matrix(dist(fit[, 1:2])))^2) / sum(matriks_jarak^2))
stress_2d
## [1] 0.2701668
#Visualisasi MDS 3D
plot_ly(dat, x = ~MDS1, y = ~MDS2, z = ~MDS3,
        type = "scatter3d", mode = "markers+text",
        text = ~Wilayah) %>%
  layout(title = "Visualisasi MDS 3D Kabupaten/Kota Provinsi Jawa Tengah berdasarkan IKK 2023",
         scene = list(xaxis = list(title = "Dimensi 1"),
                      yaxis = list(title = "Dimensi 2"),
                      zaxis = list(title = "Dimensi 3")))
stress_3d <- sqrt(sum((matriks_jarak - as.matrix(dist(fit[, 1:3])))^2) / sum(matriks_jarak^2))
stress_3d
## [1] 0.1332582
# Method 2 : K-Means Clustering (Koordinat MDS 3 Dimensi)
data_kmeans_mds <- fit[, 1:3]
colnames(data_kmeans_mds) <- c("MDS1", "MDS2", "MDS3")

# Elbow Method
fviz_nbclust(data_kmeans_mds, kmeans, method = "wss") +
  geom_vline(xintercept = 5, linetype = 2) +
  labs(subtitle = "Elbow Method - K-Means pada Koordinat MDS")

# Silhouette Method
fviz_nbclust(data_kmeans_mds, kmeans, method = "silhouette") +
  ggtitle("Silhouette Method - K-Means pada Koordinat MDS")

# k optimal = 5
set.seed(2025)
k_opt <- 5
kmeans_mds <- kmeans(data_kmeans_mds, centers = k_opt, nstart = 25)
kmeans_mds$centers
##          MDS1         MDS2       MDS3
## 1  0.71898761 -1.468218434  0.1263411
## 2  0.09762834  0.429710178 -1.1160768
## 3  1.87027041  0.850923370 -0.1085968
## 4 -2.37875976 -0.009646127 -0.0361262
## 5 -0.03983325  0.779586895  1.0269064
dat$Cluster_MDS <- as.factor(kmeans_mds$cluster)
head(dat$Cluster_MDS)
## [1] 1 1 1 1 4 1
## Levels: 1 2 3 4 5
table(dat$Cluster_MDS)
## 
##  1  2  3  4  5 
## 30 24 17 23 25
ggplot(dat, aes(x = MDS1, y = MDS2, color = Cluster_MDS, label = Wilayah)) +
  geom_point(size = 3) +
  geom_text(vjust = -1, size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = paste("Cluster Berdasarkan Koordinat MDS 3D (k =", k_opt, ")"),
       x = "Dimensi 1",
       y = "Dimensi 2",
       color = "Cluster") +
  theme_minimal() +
  theme(legend.position = "right")

fviz_cluster(kmeans_mds,
             data = data_kmeans_mds,
             geom = "point",
             main = paste("Cluster Plot pada Ruang MDS 3D (k =", k_opt, ")"))

plot_ly(dat, x = ~MDS1, y = ~MDS2, z = ~MDS3,
        type = "scatter3d", mode = "markers+text",
        text = ~Wilayah,
        color = ~Cluster_MDS) %>%
  layout(
    title = paste("Visualisasi Cluster MDS 3D (k =", k_opt, ")"),
    scene = list(
      xaxis = list(title = "Dimensi 1"),
      yaxis = list(title = "Dimensi 2"),
      zaxis = list(title = "Dimensi 3")
    )
  )
# Daftar wilayah per cluster
cluster_list_mds <- split(dat$Wilayah, dat$Cluster_MDS)
for (k in names(cluster_list_mds)) {
  cat("Cluster", k, ":\n")
  print(cluster_list_mds[[k]])
  cat("\n")
}
## Cluster 1 :
##  [1] "Kab Pandeglang"   "Kab Lebak"        "Kab Tangerang"    "Kab Serang"      
##  [5] "Kota Cilegon"     "Kota Serang"      "Kepulauan Seribu" "Bogor"           
##  [9] "Sukabumi"         "Cianjur"          "Bandung"          "Garut"           
## [13] "Kuningan"         "Cirebon"          "Sumedang"         "Indramayu"       
## [17] "Subang"           "Purwakarta"       "Karawang"         "Bekasi"          
## [21] "Bandung Barat"    "Kota Bogor"       "Kota Sukabumi"    "Kota Cirebon"    
## [25] "Kota Cimahi"      "Kota Banjar"      "Kab Cilacap"      "Kab Batang"      
## [29] "Kab Tegal"        "Kab Brebes"      
## 
## Cluster 2 :
##  [1] "Bantul"          "Kab Banyumas"    "Kab Purbalingga" "Kab Purworejo"  
##  [5] "Kab Magelang"    "Kab Boyolali"    "Kab Klaten"      "Kab Sukoharjo"  
##  [9] "Kab Wonogiri"    "Kab Karanganyar" "Kab Sragen"      "Kab Grobogan"   
## [13] "Kab Blora"       "Kab Rembang"     "Kab Pati"        "Kab Kudus"      
## [17] "Kab Jepara"      "Kab Demak"       "Kab Semarang"    "Kab Temanggung" 
## [21] "Kab Kendal"      "Kab Pekalongan"  "Kota Pekalongan" "Kota Tegal"     
## 
## Cluster 3 :
##  [1] "Kulonprogo"       "Gunungkidul"      "Tasikmalaya"      "Ciamis"          
##  [5] "Majalengka"       "Pangandaran"      "Kota Tasikmalaya" "Kab Banjarnegara"
##  [9] "Kab Kebumen"      "Kab Wonosobo"     "Kab Pemalang"     "Kab Probolinggo" 
## [13] "Kab Tuban"        "Kab Bangkalan"    "Kab Sampang"      "Kab Pamekasan"   
## [17] "Kab Sumenep"     
## 
## Cluster 4 :
##  [1] "Kota Tangerang"         "Kota Tangerang Selatan" "Sleman"                
##  [4] "Yogyakarta"             "Kota Jakarta Selatan"   "Kota Jakarta Timur"    
##  [7] "Kota Jakarta Pusat"     "Kota Jakarta Barat"     "Kota Jakarta Utara"    
## [10] "Kota Bandung"           "Kota Bekasi"            "Kota Depok"            
## [13] "Kota Magelang"          "Kota Surakarta"         "Kota Salatiga"         
## [16] "Kota Semarang"          "Kab Sidoarjo"           "Kota Blitar"           
## [19] "Kota Malang"            "Kota Mojokerto"         "Kota Madiun"           
## [22] "Kota Surabaya"          "Kota Batu"             
## 
## Cluster 5 :
##  [1] "Kab Pacitan"      "Kab Ponorogo"     "Kab Trenggalek"   "Kab Tulungagung" 
##  [5] "Kab Blitar"       "Kab Kediri"       "Kab Malang"       "Kab Lumajang"    
##  [9] "Kab Jember"       "Kab Banyuwangi"   "Kab Bondowoso"    "Kab Situbondo"   
## [13] "Kab Pasuruan"     "Kab Mojokerto"    "Kab Jombang"      "Kab Nganjuk"     
## [17] "Kab Madiun"       "Kab Magetan"      "Kab Ngawi"        "Kab Bojonegoro"  
## [21] "Kab Lamongan"     "Kab Gresik"       "Kota Kediri"      "Kota Probolinggo"
## [25] "Kota Pasuruan"
# 9. Profiling Hasil
aggregate(dat[, c("IKK", "IPM", "TPT", "PPK", "RLS")],
          by = list(Cluster = dat$Cluster_MDS),
          FUN = mean)
##   Cluster       IKK      IPM      TPT       PPK       RLS
## 1       1 1.3483333 71.06700 8.686333  9835.600  8.039333
## 2       2 1.5125000 73.70542 4.417500 15034.792  8.132083
## 3       3 2.4135294 68.91412 4.164118  9663.882  9.434118
## 4       4 0.7426087 82.20043 7.043043 16909.565 12.268696
## 5       5 1.3864000 72.31800 5.173200 11354.280 13.317600

```

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.