#Package yang digunakan

library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(corrplot)
## corrplot 0.95 loaded
library(factoextra)
## Welcome to factoextra!
## Want to learn more? See two factoextra-related books at https://www.datanovia.com/library/principal-component-methods
library(cluster)
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

1. Read csv

maternal <- read.csv("D:/STATISTIKA/SEM 5/Modern Prediksi dan Machine Learning/Tugas Project PCA, Clsuter/Maternal Health Risk Data Set.csv")
dim(maternal)
## [1] 1014    7
names(maternal)
## [1] "Age"         "SystolicBP"  "DiastolicBP" "BS"          "BodyTemp"   
## [6] "HeartRate"   "RiskLevel"
head(maternal)
##   Age SystolicBP DiastolicBP    BS BodyTemp HeartRate RiskLevel
## 1  25        130          80 15.00       98        86 high risk
## 2  35        140          90 13.00       98        70 high risk
## 3  29         90          70  8.00      100        80 high risk
## 4  30        140          85  7.00       98        70 high risk
## 5  35        120          60  6.10       98        76  low risk
## 6  23        140          80  7.01       98        70 high risk
str(maternal)
## 'data.frame':    1014 obs. of  7 variables:
##  $ Age        : int  25 35 29 30 35 23 23 35 32 42 ...
##  $ SystolicBP : int  130 140 90 140 120 140 130 85 120 130 ...
##  $ DiastolicBP: int  80 90 70 85 60 80 70 60 90 80 ...
##  $ BS         : num  15 13 8 7 6.1 7.01 7.01 11 6.9 18 ...
##  $ BodyTemp   : num  98 98 100 98 98 98 98 102 98 98 ...
##  $ HeartRate  : int  86 70 80 70 76 70 78 86 70 70 ...
##  $ RiskLevel  : chr  "high risk" "high risk" "high risk" "high risk" ...
colSums(is.na(maternal))
##         Age  SystolicBP DiastolicBP          BS    BodyTemp   HeartRate 
##           0           0           0           0           0           0 
##   RiskLevel 
##           0

2a. EKSPLORASI DATA

# Memilih variabel numerik untuk analisis
data.pca = maternal[, c("Age", "SystolicBP","DiastolicBP",
                        "BS", "BodyTemp", "HeartRate")]
#Statistik deskriptif
str(data.pca)
## 'data.frame':    1014 obs. of  6 variables:
##  $ Age        : int  25 35 29 30 35 23 23 35 32 42 ...
##  $ SystolicBP : int  130 140 90 140 120 140 130 85 120 130 ...
##  $ DiastolicBP: int  80 90 70 85 60 80 70 60 90 80 ...
##  $ BS         : num  15 13 8 7 6.1 7.01 7.01 11 6.9 18 ...
##  $ BodyTemp   : num  98 98 100 98 98 98 98 102 98 98 ...
##  $ HeartRate  : int  86 70 80 70 76 70 78 86 70 70 ...
summary(data.pca)
##       Age          SystolicBP     DiastolicBP           BS        
##  Min.   :10.00   Min.   : 70.0   Min.   : 49.00   Min.   : 6.000  
##  1st Qu.:19.00   1st Qu.:100.0   1st Qu.: 65.00   1st Qu.: 6.900  
##  Median :26.00   Median :120.0   Median : 80.00   Median : 7.500  
##  Mean   :29.87   Mean   :113.2   Mean   : 76.46   Mean   : 8.726  
##  3rd Qu.:39.00   3rd Qu.:120.0   3rd Qu.: 90.00   3rd Qu.: 8.000  
##  Max.   :70.00   Max.   :160.0   Max.   :100.00   Max.   :19.000  
##     BodyTemp        HeartRate   
##  Min.   : 98.00   Min.   : 7.0  
##  1st Qu.: 98.00   1st Qu.:70.0  
##  Median : 98.00   Median :76.0  
##  Mean   : 98.67   Mean   :74.3  
##  3rd Qu.: 98.00   3rd Qu.:80.0  
##  Max.   :103.00   Max.   :90.0
# Visualisasi statistik deskriptif
par(mfrow = c(2,3))
hist(data.pca$Age, main = "Distribusi Age", xlab = "Age")
hist(data.pca$SystolicBP, main = "Distribusi Systolic BP", xlab = "Systolic BP")
hist(data.pca$DiastolicBP, main = "Distribusi Diastolic BP", xlab = "Diastolic BP")
hist(data.pca$BS, main = "Distribusi Blood Sugar", xlab = "BS")
hist(data.pca$BodyTemp, main = "Distribusi Body Temperature", xlab = "BodyTemp")
hist(data.pca$HeartRate, main = "Distribusi Heart Rate", xlab = "HeartRate")

par(mfrow = c(1,1))

# Boxplot
boxplot(data.pca, main = "Boxplot Variabel Kesehatan",las = 2)

# Matriks korelasi
kor = cor(data.pca)

# Visualisasi korelasi
corrplot(kor, method = "color", type = "upper",tl.col = "black")

corrplot(kor)

2b. PREPROCESSING

# Standardisasi data
maternalscale = scale(data.pca)
head(maternalscale)
##               Age SystolicBP DiastolicBP         BS   BodyTemp  HeartRate
## [1,] -0.361559706  0.9129458   0.2548970  1.9049502 -0.4849762  1.4462425
## [2,]  0.380589164  1.4563085   0.9750574  1.2976993 -0.4849762 -0.5318251
## [3,] -0.064700158 -1.2605050  -0.4652634 -0.2204279  0.9734042  0.7044671
## [4,]  0.009514729  1.4563085   0.6149772 -0.5240533 -0.4849762 -0.5318251
## [5,]  0.380589164  0.3695831  -1.1854238 -0.7973162 -0.4849762  0.2099502
## [6,] -0.509989480  1.4563085   0.2548970 -0.5210171 -0.4849762 -0.5318251
apply(maternalscale, 2, mean)
##           Age    SystolicBP   DiastolicBP            BS      BodyTemp 
## -3.101660e-17 -6.714943e-17  4.532743e-16 -2.549239e-16  2.120426e-15 
##     HeartRate 
##  1.743124e-16
apply(maternalscale, 2, sd)
##         Age  SystolicBP DiastolicBP          BS    BodyTemp   HeartRate 
##           1           1           1           1           1           1
# Boxplot sebelum standardisasi
boxplot(data.pca, main = "Boxplot Data Sebelum Standardisasi", las = 2)

# Boxplot setelah standardisasi
boxplot(maternalscale, main = "Boxplot Data Setelah Standardisasi", las = 2)

# 2c. PRINCIPAL COMPONENT ANALYSIS

# Perhitungan nilai eigen
cbind( eigen(cov(maternalscale))$values, eigen(cor(data.pca))$values)
##           [,1]      [,2]
## [1,] 2.6078934 2.6078934
## [2,] 1.1443812 1.1443812
## [3,] 0.8370499 0.8370499
## [4,] 0.7063345 0.7063345
## [5,] 0.4925435 0.4925435
## [6,] 0.2117975 0.2117975
# Analisis PCA
hasilpca = prcomp( data.pca, scale. = TRUE)
summary(hasilpca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5    PC6
## Standard deviation     1.6149 1.0698 0.9149 0.8404 0.70181 0.4602
## Proportion of Variance 0.4346 0.1907 0.1395 0.1177 0.08209 0.0353
## Cumulative Proportion  0.4346 0.6254 0.7649 0.8826 0.96470 1.0000
# Nilai eigen
eigenvalue = hasilpca$sdev^2
eigenvalue
## [1] 2.6078934 1.1443812 0.8370499 0.7063345 0.4925435 0.2117975
# Proporsi varians
proporsi = eigenvalue / sum(eigenvalue)
proporsi
## [1] 0.43464890 0.19073020 0.13950832 0.11772242 0.08209058 0.03529959
# Persentase varians
persen = proporsi * 100
persen
## [1] 43.464890 19.073020 13.950832 11.772242  8.209058  3.529959
# Persentase kumulatif
kumulatif = cumsum(persen)
kumulatif
## [1]  43.46489  62.53791  76.48874  88.26098  96.47004 100.00000
# Scree Plot
plot(hasilpca$sdev^2, type = "o", main = "Scree Plot",
     xlab = "Komponen Utama", ylab = "Eigenvalue")

# Tabel perbandingan PCA
Perbandingan = data.frame(
  Lambda = round(hasilpca$sdev^2, 4),
  Pctn = round(hasilpca$sdev^2 /sum(hasilpca$sdev^2) * 100, 3),
  CumPctn = cumsum(round(hasilpca$sdev^2 / sum(hasilpca$sdev^2) * 100, 2)
  )
)
Perbandingan
##   Lambda   Pctn CumPctn
## 1 2.6079 43.465   43.46
## 2 1.1444 19.073   62.53
## 3 0.8370 13.951   76.48
## 4 0.7063 11.772   88.25
## 5 0.4925  8.209   96.46
## 6 0.2118  3.530   99.99
# Loading PC1-PC3
hasilpca$rotation[,1:3]
##                    PC1        PC2        PC3
## Age          0.4360752 -0.1729019  0.2404386
## SystolicBP   0.5296016  0.1128996 -0.2415842
## DiastolicBP  0.5225675  0.1230416 -0.2990961
## BS           0.4255290 -0.3528904 -0.1153507
## BodyTemp    -0.2735090 -0.4293197 -0.8093439
## HeartRate    0.0200401 -0.7958469  0.3549994
# Skor komponen utama
maternalpca = data.frame(maternalscale %*% hasilpca$rotation[,1:3])

# Visualisasi PCA 3D
plot_ly(maternalpca, x = ~PC1, y = ~PC2,z = ~PC3)
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

2d. K-MEANS CLUSTERING

# Menentukan jumlah cluster dengan WSS
wssplot = function(data, nc=10, seed=1234){
  wss = (nrow(data)-1) *
    sum(apply(data,2,var))
  
  for (i in 2:nc){
    set.seed(seed)
    wss[i] <- sum(
      kmeans(
        data,
        centers=i
      )$withinss
    )
  }
  
  plot(
    1:nc,
    wss,
    type="b",
    xlab="Number of Clusters",
    ylab="Within groups sum of squares",
    main="WSS Plot"
  )
}
wssplot(maternalscale)

# Menentukan jumlah cluster dengan Silhouette
fviz_nbclust(maternalscale, kmeans, method="silhouette")

# Berdasarkan WSS dan Silhouette,
# jumlah cluster optimal = 4

# K-Means final
set.seed(123)
clustering = kmeans(maternalscale, centers=4)
clustering
## K-means clustering with 4 clusters of sizes 177, 247, 181, 409
## 
## Cluster means:
##          Age SystolicBP DiastolicBP         BS   BodyTemp   HeartRate
## 1 -0.5858815 -0.6922766  -0.6288252 -0.2305487  2.0321718  0.28887736
## 2 -0.2828378 -1.0046626  -1.1093258 -0.4743792 -0.4530926  0.01574723
## 3  1.1731714  0.9786897   0.9810256  1.8806266 -0.4181002  0.39915075
## 4 -0.0948216  0.4732073   0.5079216 -0.4460015 -0.4207932 -0.31116661
## 
## Clustering vector:
##    [1] 3 3 1 4 2 4 4 1 4 3 2 4 4 4 3 4 3 4 1 2 3 4 2 2 4 2 4 2 4 2 4 4 4 4 2 1 4
##   [38] 2 4 2 4 4 2 4 4 4 4 2 2 4 4 2 2 4 2 4 4 4 2 2 4 2 2 2 2 4 1 4 1 1 1 4 4 3
##   [75] 3 4 2 3 4 2 2 2 2 4 2 2 4 2 4 4 4 3 1 4 1 4 4 4 4 2 4 4 3 4 4 1 3 4 3 4 3
##  [112] 1 2 3 3 4 1 3 4 3 3 3 4 3 3 4 3 3 4 3 3 3 4 3 3 1 3 3 3 1 4 1 1 1 1 1 2 2
##  [149] 2 3 2 4 1 4 2 2 4 3 3 4 4 2 2 4 2 4 3 3 2 4 2 1 4 4 4 4 4 3 3 2 1 3 3 3 4
##  [186] 2 4 2 4 1 4 1 1 3 2 4 2 4 4 1 4 4 2 4 2 4 3 4 4 2 3 2 4 2 1 4 2 3 2 4 4 1
##  [223] 4 4 2 4 2 4 3 4 3 4 4 1 3 1 4 1 1 1 1 1 2 1 2 3 2 4 1 3 2 2 4 4 4 4 4 2 2
##  [260] 4 2 4 3 3 2 4 2 1 4 4 4 4 4 3 3 2 1 3 3 3 4 2 4 2 4 1 4 1 1 4 2 4 2 4 4 1
##  [297] 4 4 2 4 2 4 3 4 4 2 4 2 4 2 4 2 2 4 2 4 2 3 3 3 3 4 2 4 4 2 2 3 2 3 2 3 4
##  [334] 2 3 4 2 1 1 1 1 1 1 2 2 3 2 4 1 4 2 2 4 4 4 4 4 2 2 4 2 4 3 3 4 4 4 4 3 3
##  [371] 2 1 3 3 3 4 2 4 2 4 1 4 1 1 4 2 4 2 4 4 1 4 4 2 4 2 4 3 4 4 2 4 2 4 2 4 4
##  [408] 2 4 2 4 2 4 3 1 3 1 3 4 1 4 2 4 2 4 2 3 4 4 2 4 2 3 4 4 3 3 2 4 2 4 4 1 4
##  [445] 3 2 4 2 4 3 4 3 4 2 4 2 4 4 1 4 3 2 4 2 4 3 4 3 4 1 4 3 1 4 1 4 2 4 3 4 4
##  [482] 2 3 2 4 2 1 4 2 3 2 4 4 1 4 4 2 4 2 4 3 4 3 4 1 2 3 1 4 1 1 1 1 1 2 2 1 4
##  [519] 2 3 2 4 4 1 4 4 2 4 2 4 3 4 3 4 1 2 2 4 3 3 2 4 2 1 4 4 4 4 4 3 3 2 1 3 3
##  [556] 3 4 2 4 2 4 1 4 1 1 4 4 4 2 4 4 1 4 4 2 4 2 4 3 4 4 2 4 3 3 2 4 2 1 4 4 4
##  [593] 4 4 3 3 2 1 3 3 3 4 3 3 3 4 2 4 2 4 1 4 1 1 3 2 4 2 4 4 1 4 4 2 4 2 4 3 4
##  [630] 4 2 3 2 4 2 1 3 3 1 2 3 1 4 2 3 2 4 4 1 4 4 2 4 2 4 3 4 3 3 4 1 2 4 2 2 2
##  [667] 2 2 2 2 2 2 2 2 2 3 1 4 4 1 3 3 4 3 4 1 3 1 3 3 4 4 4 2 2 4 2 4 2 4 4 2 2
##  [704] 2 2 2 4 4 4 4 3 1 1 1 1 1 1 1 2 4 4 1 4 2 2 4 1 4 1 4 4 3 1 1 1 1 1 1 1 2
##  [741] 4 4 1 4 2 2 4 1 4 1 1 1 4 4 3 3 4 2 3 4 2 2 2 2 4 2 2 4 2 4 4 4 3 1 4 1 4
##  [778] 4 4 4 2 4 4 4 2 4 2 3 1 1 4 4 4 2 4 2 4 3 4 4 4 4 3 4 4 4 4 2 4 4 3 1 1 4
##  [815] 4 4 1 3 3 3 1 2 4 4 2 4 2 1 4 4 4 4 4 3 1 4 4 4 2 1 4 2 4 1 4 4 4 1 2 2 2
##  [852] 2 4 4 4 4 3 1 1 1 1 1 1 1 2 4 4 1 4 2 2 4 1 4 1 1 4 4 4 2 2 4 2 2 4 2 2 4
##  [889] 3 4 1 1 1 2 4 3 4 2 2 4 2 4 1 4 4 2 4 2 4 4 2 1 1 1 1 1 2 2 4 2 2 4 1 4 4
##  [926] 2 4 2 2 2 4 2 4 2 1 4 4 4 4 2 1 4 2 4 2 4 1 4 1 4 4 1 4 4 4 3 4 1 3 1 3 3
##  [963] 3 1 4 3 3 3 3 3 1 1 3 1 4 3 3 3 4 3 3 4 4 3 1 4 1 1 1 1 1 3 1 3 3 3 1 4 3
## [1000] 3 1 3 3 3 1 1 1 4 3 2 3 3 3 1
## 
## Within cluster sum of squares by cluster:
## [1] 542.0906 581.4040 650.1096 962.8363
##  (between_SS / total_SS =  55.0 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# Jumlah anggota setiap cluster
table(clustering$cluster)
## 
##   1   2   3   4 
## 177 247 181 409
# Centroid setiap cluster
clustering$centers
##          Age SystolicBP DiastolicBP         BS   BodyTemp   HeartRate
## 1 -0.5858815 -0.6922766  -0.6288252 -0.2305487  2.0321718  0.28887736
## 2 -0.2828378 -1.0046626  -1.1093258 -0.4743792 -0.4530926  0.01574723
## 3  1.1731714  0.9786897   0.9810256  1.8806266 -0.4181002  0.39915075
## 4 -0.0948216  0.4732073   0.5079216 -0.4460015 -0.4207932 -0.31116661
# Nilai silhouette rata-rata
sil = silhouette(clustering$cluster, dist(maternalscale))
mean(sil[,3])
## [1] 0.3063386
plot(sil)

# VISUALISASI HASIL CLUSTER PADA RUANG PCA
pca = prcomp(maternalscale)
maternalpca = data.frame(maternalscale %*% pca$rotation)
maternalpca$cluster = factor(clustering$cluster)
plot_ly(x=~PC1, y=~PC2, z=~PC3, data=maternalpca, color=~cluster)
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
# PROFIL CLUSTER
maternal2 = cbind(data.pca, cluster=clustering$cluster)

# Rata-rata setiap variabel berdasarkan cluster
aggregate(data.pca, by=list(Cluster=clustering$cluster),FUN=mean)
##   Cluster      Age SystolicBP DiastolicBP        BS  BodyTemp HeartRate
## 1       1 21.97740   100.4576    67.72881  7.966667 101.45198  76.63842
## 2       2 26.06073    94.7085    61.05668  7.163603  98.04372  74.42915
## 3       3 45.67956   131.2099    90.08287 14.919890  98.09171  77.53039
## 4       4 28.59413   121.9071    83.51345  7.257066  98.08802  71.78484
# Boxplot profil setiap cluster
par(mfrow=c(2,3))
boxplot(Age ~ cluster, maternal2, main="Age")
boxplot(SystolicBP ~ cluster, maternal2, main="Systolic BP")
boxplot(DiastolicBP ~ cluster, maternal2,main="Diastolic BP")
boxplot(BS ~ cluster, maternal2, main="Blood Sugar")
boxplot(BodyTemp ~ cluster, maternal2, main="Body Temperature")
boxplot(HeartRate ~ cluster,  maternal2, main="Heart Rate")

par(mfrow=c(1,1))

TAMBAHAN: CLUSTER DAN RISK LEVEL

maternal2 = cbind(maternal, cluster=clustering$cluster)
# Tabel frekuensi
table(maternal2$cluster, maternal2$RiskLevel)
##    
##     high risk low risk mid risk
##   1        62       37       78
##   2        10      168       69
##   3       144        5       32
##   4        56      196      157
# Proporsi RiskLevel dalam setiap cluster
prop.table(table(maternal2$cluster, maternal2$RiskLevel), margin=1)
##    
##      high risk   low risk   mid risk
##   1 0.35028249 0.20903955 0.44067797
##   2 0.04048583 0.68016194 0.27935223
##   3 0.79558011 0.02762431 0.17679558
##   4 0.13691932 0.47921760 0.38386308