07-UL.Rmd

# Unsupervised Learning

```{r setup, include=FALSE}
# clear-up the environment
rm(list = ls())

# chunk options
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  comment = "#>"
)

# scientific notation
options(scipen = 9999)
```

```{r echo=FALSE, warning=FALSE, message=FALSE}
# libraries
library(tidyverse)
library(FactoMineR)
library(factoextra)
library(rsample)
library(recipes)
library(plotly)
```

## PCA: Dimensionality Reduction

### **Permasalahan apa yang terdapat pada data berdimensi tinggi?**

  - menyulitkan pengolahan data
  - memerlukan komputasi yang besar
  - tidak efisien secara waktu

###  **Perbedaan membuat PCA dengan menggunakan fungsi `prcomp()` dan `PCA()` dari library `FactoMineR`?**

Fungsi untuk membuat biplot di R:

- `biplot(prcomp())` -> base R
- `plot.PCA(PCA())` -> package `FactoMineR`

Kelebihan ketika membuat PCA dengan menggunakan fungsi `PCA()` dari library `FactoMineR` adalah bisa membuat biplot lebih spesifik:
  
- memisahkan dua grafik yang terdapat pada biplot yaitu **individual factor map** dan **variables factor map**
- mengkombinasikan antara variabel **numerik dan kategorik** dengan menggunakan fungsi `plot.PCA()`. 

### **Apakah terdapat best practice dalam menentukan jumlah PC yang digunakan pada PCA? **

   Penentuan jumlah PC yang digunakan bergantung pada kebutuhan analisa yang dilakukan. Namun, kembali pada tujuan awal melakukan PCA, yaitu untuk mereduksi dimensi supaya analisis lanjutan yang dilakukan memiliki waktu yang relatif cepat dan ruang penyimpanan yang lebih efisien. Sehingga, seringkali seorang analis menentapkan threshold lebih dari 70-75% informasi. Maksudnya jumlah PC yang digunakan adalah jumlah PC yang sudah merangkum kurang lebih 70-75% informasi. Namun, threshold tersebut sifatnya tidak mutlak artinya disesuaikan dengan kebutuhan analisis dan bisnis. Berikut merupakan link eksternal yang dapat dijadikan sebagai bahan referensi [How many components can I retrieve in principal component analysis?](https://www.researchgate.net/post/How_many_components_can_I_retrieve_in_principal_component_analysis).

### **Bagaimana implementasi PCA pada data pre-processing?**

Berikut ini adalah implementasi PCA pada tahapan data pre-processing dengan menggunakan data `attrition`.

```{r}
attrition <- read.csv("data/06-UL/attrition.csv")
str(attrition)
```

Sebelum melakukan PCA terlebih dahulu dilakukan cross validation, yaitu membagi data menjadi **training set** untuk proses pemodelan dan **testing set** untuk melakukan evaluasi. Namun, data train dan data test tidak langsung dimasukkan ke dalam sebuah objek melainkan dilakukan PCA terlebih dahulu.

Cross validation akan dilakukan dengan menggunakan fungsi `initial_split()` dari library `rsample`. Fungsi tersebut akan melakukan proses sampling untuk cross validation dengan metode **stratified random sampling**, sehingga proporsi target variabel pada data awal, akan dipertahankan baik pada training set maupun testing set.

```{r}
set.seed(417)
splitted <- initial_split(data = attrition,
                          prop = 0.8,
                          strata = "attrition")
splitted
```

Melakukan tahapan data preparation yang didalamnya termasuk melakukan PCA. Data preparation yang akan dilakukan adalah menghapus variabel yang dianggap tidak berpengaruh, membuang variabel yang variansinya mendekati 0 (tidak informatif), melakukan scaling, dan melakukan PCA. Proses yang dilakukan pada tahapan data preparation akan dilakukan dengan menggunakan fungsi dari library `recipes`, yaitu:

- `step_rm()` untuk menghapus variabel
- `step_nzv()` untuk membuang variabel yang variansinya mendekati 0
- `step_center()` dan `step_scale()` untuk melakukan scaling
- `step_pca()` untuk melakukan PCA

```{r}
rec <- recipe(attrition ~ ., training(splitted)) %>% 
  step_rm(employee_count, employee_number) %>%
  step_nzv(all_predictors()) %>% 
  step_center(all_numeric()) %>% 
  step_scale(all_numeric()) %>% 
  step_pca(all_numeric(), threshold = 0.8) %>% 
  prep()
```

Setelah mendefinisikan proses data preparation pada objek `rec`, selanjutnya proses tersebut diterapkan ke data train menggunakan fungsi `juice()` dan ke data test menggunakan fungsi `bake()` dari library `recipes`.

```{r}
train <- juice(rec)
head(train, 5)
```

```{r}
test <- bake(rec, testing(splitted))
head(test, 5)
```

Dari output di atas diketahui bahwa variabel numerik sudah berbentuk sebuah PC. Selanjutnya, data sudah siap dilanjutkan ke tahap modeling.

### **Bagaimana penerapan PCA di industri?**

  PCA pada industri lebih sering digunakan untuk data preparation sama halnya seperti scaling, feature engineering, ataupun feature selection. PCA digunakan untuk mereduksi data berdimensi besar besar menjadi lebih kecil, secara sederhana dapat dikatakan mengurangi jumlah kolom pada data. Walaupun begitu, PCA tetap mempertahankan informasi dari semua variabel. Sebelum mereduksi dimensi, PCA akan merangkum terlebih dahulu semua informasi yang terdapat pada setiap variabel ke dalam bentuk PC, PC tersebut yang nantinya akan direduksi (dikurangi) dimensinya. Oleh karena itu, variabel yang digunakan jumlahnya tetap sama seperti data awal, hanya informasi (variansinya) saja yang berkurang. Berikut merupakan link eksternal yang dapat dijadikan sebagai bahan referensi [An Application of PCA](https://web.cs.ucdavis.edu/~vemuri/papers/pcaVisualization.pdf).

   Contoh permasalahan yang sering ditemui adalah klasifikasi dokumen. Saat ini semua administrasi dilakukan secara online/elektronik (tidak manual), adakalanya seorang nasabah/pelamar/customer harus melakukan upload dokumen. Sebelum adanya klasifikasi dokumen, pemeriksaan kebenaran dokumen dilakukan secara manual sehingga membutuhkan waktu yang cukup lama dan kapasitas penyimpanan yang relatif besar karena aplikasi tidak mampu memilah mana dokumen yang sudah sesuai dan mana yang belum. Namun, permasalahan tersebut sudah mampu terjawab dengan adanya klasifikasi dokumen. Data untuk klasifikasi dokumen adalah data image yang jika dilakukan proses klasifikasi akan memerlukan komputasi yang relatif lama dibandingkan data tabular biasa. Oleh karena itu, perlu dilakukan PCA untuk mereduksi dimensi data image tersebut supaya komputasi saat proses klasifikasi bisa menjadi lebih cepat. Berikut merupakan link eksternal yang dapat dijadikan sebagai bahan referensi [Image Compression with PCA in R](https://www.r-bloggers.com/image-compression-with-pca-in-r/).

## PCA: Visualization

### **Apakah biplot dapat menampilkan PC lain selain PC1 dan PC2?**

  Bisa, tetapi informasi yang dijelaskan menjadi berkurang, karena secara default PC1 dan PC2 merangkum informasi paling banyak. Berikut contoh membuat biplot dengan menggunakan PC lain (selain PC1 dan PC2):
  
```{r}
head(USArrests)
```

Membuat PCA dari data `USArrests` dengan menggunakan fungsi `prcomp()`.

```{r}
pca_us <- prcomp(USArrests, scale = T)
```

Membuat visualisasi dari hasil PCA dengan menggunakan fungsi `biplot()`.

```{r}
# parameter `choices` dapat diganti sesuai PC yang ingin dibuat, secara default menggunakan  PC1 dan PC2 (choices = 1:2)
biplot(pca_us, choices = 2:3)
```
   
### **Apakah kita dapat memvisualisasikan biplot dengan 3 dimensi?**

Untuk menampilkan biplot dengan 3 dimensi dapat menggunakan function `plot_ly()` dari package `plotly`. Berikut ini akan dicontohkan memvisualisasikan biplot dari PC1, PC2, PC3 dan juga akan dibedakan setiap titik observasi dengan cluster nya. Sebelum masuk ke visualisasi, akan dicari terlebih dahulu cluster untuk setiap observasi.

```{r}
# Read data in
whiskies <- read.csv("data/06-UL/whiskies.txt")

# Distillery column is the name of each whisky
rownames(whiskies) <- whiskies[,"Distillery"]

# remove RowID, Postcode, Latitude and Longitude
whiskies <- whiskies[,3:14]

# k-means clustering
whi_km <- kmeans(scale(whiskies), 4)
```

Setelah menggunakan `kmeans()` untuk mendapatkan cluster, berikutnya kita lakukan PCA dan membentuk PC yang diperoleh dalam bentuk data frame.

```{r}
whis.pca <- PCA(whiskies, graph = F, scale.unit = T)
df_pca <- data.frame(whis.pca$ind$coord) %>% 
          bind_cols(cluster = as.factor(whi_km$cluster))
head(df_pca)
```
 
Visualisasikan PC dan membedakan warna tiap observasi berdasarkan clusternya.

```{r eval=FALSE}
plot_ly(df_pca, x = ~Dim.1, y = ~Dim.2, z = ~Dim.3, color = ~cluster)
```

### **Bagaimana implementasi visualisasi PCA menggunakan package `factoextra`?**

Kita akan mencoba melakukan visualisasi pada data setelah dilakukan PCA dengan menggunakan data `loan`.

```{r}
loan <- read.csv("data/06-UL/loan2017Q4.csv")
str(loan)
```

Sebelum melakukan PCA kita akan melakukan tahapan data preparation terlebih dahulu dengan membuang variabel `initial_list_status`, `home_ownership`, dan `not_paid` karena visualisasi yang akan dibuat tidak memerlukan insight dari ketiga variabel tersebut.

```{r}
loan_clean <- loan %>% 
  select(-c(initial_list_status, home_ownership, not_paid))
```

Membuat PCA dengan menggunakan fungsi `PCA()` dari library `FactoMineR`. Parameter yang digunakan adalah:

- `ncp`: Jumlah PC yang akan dihasilkan. Secara default fungsi `PCA()` hanya akan menampilkan 5 PC awal (5 PC yang merangkum informasi paling banyak)
- `quali.sup`: Nomor kolom dari variabel kategorik
- `graph`: Sebuah logical value. `T` akan menampilkan hasil visualisasi, `F` tidak menampilkan hasil visualisasi. Secara default fungsi `PCA()` akan langsung menampilkan hasil visualisasi

```{r}
pca_loan <- PCA(loan_clean,
                ncp = 10,
                quali.sup = c(1, 6, 7),
                graph = F)
```

Setelah membuat PCA, selanjutnya adalah membuat visualisasi dari hasil PCA. Kita akan membuat individual plot menggunakan fungsi `fviz_pca_ind()` dari library `factoextra`. Parameter yang digunakan adalah:

- Objek hasil PCA
- `habillage`: Nomor kolom dari variabel kategorik, setiap individu akan dibedakan berdasarkan variabel kategori yang dipilih
- `select.ind`: Jumlah individu dengan kontribusi tertinggi yang ingin dilihat

```{r}
fviz_pca_ind(pca_loan,
             habillage = 6,
             select.ind = list(contrib = 10))
```

Plot individu di atas hanya menampilkan 10 observasi yang memberikan informasi tertinggi terhadap PC1 dan PC2. Namun, terdapat lebih dari 10 titik observasi yang terdapat pada plot di atas karena terdapat titik observasi yang merupakan titik pusat dari tiap status verifikasi.

### **Bagaimana cara PCA menentukan outlier dari individual plot dan langsung ekstrak daftar outlier menjadi sebuah vector?**

Misalkan kita menggunakan individual plot dari `pca_loan` untuk melihat 5 observasi yang tergolong outlier, sebagai berikut:

```{r}
plot.PCA(pca_loan, choix = "ind", select = "contrib 5")
```

Tujuan kita adalah untuk mengekstrak label outlier yang ditampilkan pada individual plot di atas. Berdasarkan dokumentasi fungsi `plot.PCA()`, outlier ditentukan melalui jumlah kuadrat koordinat (Sum of Squares `SS`) tertinggi dari dimensi yang dipilih pada individual plot. Secara default, PC yang digunakan adalah PC1 (`Dim.1`) dan PC2 (`Dim.2`).

```{r}
coord_mat <- pca_loan$ind$coord
data.frame("outlier_name" = row.names(coord_mat),
           coord_mat) %>% 
  mutate(SS = Dim.1^2 + Dim.2^2) %>% # sesuaikan dengan PC yang dipilih
  arrange(-SS) %>% 
  head(5) %>% # sesuaikan dengan jumlah outlier yang diambil
  pull(outlier_name)
```

## Clustering

### **Bagaimana best practice dalam penentuan jumlah cluster?**

   Fungsi `kmeans()` tidak dapat menentukan jumlah cluster secara otomatis. Jumlah cluster tetap ditentukan oleh user berdasarkan kebutuhan bisnis. Namun, secara statistik penentuan jumlah cluster dapat dilakukan berdasarkan penurunan Within Sum of Square (WSS). Secara sederhana, penurunan WSS dapat divisualisasikan dengan menggunakan fungsi `fviz_nbclust()` dari library `factoextra`. Berikut contoh memvisualisasikan penurunan WSS dengan menggunakan data `USArrests`:
   
```{r}
head(USArrests, 6)

# scaling data
USArrests_scale <- scale(USArrests)
```

Melakukan visualisasi penurunan WSS dengan menggunakan fungsi `fviz_nbclust()` dari library `factoextra`.

```{r}
set.seed(100)
fviz_nbclust(USArrests_scale, method = "wss", kmeans)
```

Jumlah cluster yang dipilih adalah jumlah cluster yang ketika dilakukan penambahan cluster sudah tidak mengakibatkan penurunan WSS yang signifikan (pada grafik bentuknya landai), kemudian disesuaikan dengan kebutuhan bisnis pada industri.

### **Bagaimana implementasi visualisasi K-means clustering menggunakan package `factoextra`?**

Kita akan mencoba melakukan visualisasi hasil clustering dengan menggunakan data `USArrests`.

```{r}
head(USArrests, 5)

# scaling data `USArrests`
USArrests_scale <- scale(USArrests)
```

Menentukan jumlah cluster yang akan dibuat berdasarkan penurunan WSS, dengan menggunakan fungsi `fviz_nbclust()` dari library `factoextra`.

```{r}
set.seed(100)
fviz_nbclust(USArrests_scale, method = "wss", kmeans)
```

Melakukan k-means clustering dengan jumlah cluster 5 berdasarkan hasil penurunan wss di atas menggunakan fungsi `kmeans()`. 

```{r}
set.seed(100)
USArrests_cl <- kmeans(USArrests_scale, centers = 5)
```

Membuat visualisasi hasil cluster dengan menggunakan fungsi `fviz_cluster()` dari library `factoextra`.

```{r}
# `USArrests` bisa diganti dengan yang sudah dilakukan scaling `USArrests_scale` 
fviz_cluster(USArrests_cl, USArrests)
```

Mengkombinasikan visualisasi hasil clustering dengan PCA. Untuk melakukan hal tersebut kita harus menambahkan kolom cluster pada data `USArrests`.

```{r}
USArrests_clean <- USArrests %>% 
  mutate(cluster = USArrests_cl$cluster)

head(USArrests_clean, 5)
```

Mengubah nama baris yang awalnya berupa indeks menjadi nama negara sesuai dengan data `USArrests`.

```{r}
rownames(USArrests_clean) <- rownames(USArrests)
head(USArrests_clean, 5)
```

Membuat PCA terlebih dahulu untuk mengkombinasikan visualisasi hasil clustering dengan PCA dengan menggunakan `PCA()`.

```{r}
pca_USArrest <- PCA(USArrests_clean, quali.sup = 5, graph = F)
```

Mengkombinasikan visualisasi hasil clustering dan PCA menggunakan fungsi `fviz_pca_biplot()` dari library `factoextra`. Parameter yang digunakan adalah:

- Objek hasil PCA
- `habillage`: Nomor kolom dari variabel kategorik, setiap individu akan dibedakan berdasarkan variabel kategori yang dipilih.
- `addEllipses`: Sebuah logical value. `T` akan menambah elips untuk ssetiap cluster, `F` sebaliknya. Secara default fungsi `fviz_pca_biplot()` tidak akan menambah elips pada plot individu.
  
```{r}
fviz_pca_biplot(pca_USArrest, habillage = 5, addEllipses = T)
```

Dari plot di atas terlihat bahwa antar cluster saling tumpang tindih, namun kenyataannya antar cluster pasti memiliki observasi/individu yang unik. Hal tersebut terjadi karena kita mencoba untuk memvisualisasikan cluster yang dibentuk dari 4 dimensi menjadi 2 dimensi saja.

### **Bagaimana memberikan penamaan label dari hasil clustering?**

Algoritma clustering hanya memberikan label 1 sampai $k$, di mana $k$ adalah banyaknya cluster yang terbentuk. Untuk melakukan cluster profiling atau pemberian nama cluster harus berdasarkan hasil interpretasi manusia. Alat bantu cluster profiling yang dapat digunakan adalah:

- Melalui visualisasi PCA - biplot menggunakan `fviz_pca_biplot()`
- Melihat rangkuman nilai rata-rata (centroid) untuk masing-masing cluster

Sebagai contoh, dari hasil sebelumnya, didapatkan 5 clusters untuk mengelompokkan negara bagian US berdasarkan tingkat kriminalitas dan besar populasinya. Berikut interpretasi clusternya:

```{r}
USArrests_clean %>% 
  group_by(cluster) %>% 
  summarize_all(mean)
```

- Cluster 1: HIGH crime rate, but MODERATE population -> "HIGH-MODERATE"
- Cluster 2: LOW crime rate, but HIGH population -> "LOW-HIGH"
- Cluster 3: LOWEST crime rate and population -> "LOWEST"
- Cluster 4: HIGHEST crime rate and population -> "HIGHEST"
- Cluster 5: MODERATE crime rate and population -> "MODERATE"

Menyiapkan dataframe berisi pemetaan dari nomor cluster menjadi label yang diinginkan:

```{r}
mapping_label <- data.frame(cluster = 1:5,
                            label = c("HIGH-MODERATE", "LOW-HIGH", 
                                      "LOWEST", "HIGHEST",
                                      "MODERATE"))

mapping_label
```

Melakukan join table `USArrests_clean` dan `mapping_label` dengan fungsi `join()` dari library `plyr`:

```{r}
USArrests_clean %>% 
  rownames_to_column("State") %>% 
  plyr::join(mapping_label) %>% 
  select(State, label)
```

## Mathematics Concept

### Principal Component Analysis

**1. Tujuan: membuat axis baru yang dapat merangkum informasi data**

```{r echo=FALSE, out.width="60%"}
knitr::include_graphics("assets/06-UL/pca1.png")
```

Kita akan menggunakan data dummy yang disimpan pada objek A, yang terdiri dari variabel x dan y.

```{r}
set.seed(100)
x <- runif(200)
A <- data.frame(x=x, y=-x+runif(100, 1.05, 1.25))
head(A)
```

Sebelum menghitung matriks variance-covariance harus dilakukan scaling terlebih dahulu, karena nilai variance dan covariance sangat bergantung pada interval nilai dari data.

```{r}
A <- scale(A, center = T)
head(A)
```

Membuat visualisasi data yang disimpan pada objek A.

```{r}
plot(A, cex=0.4)
```

Variabel x dan y berkorelasi secara negatif, karena variabel y dibentuk dari operasi matematika dari variabel x (dengan kata lain variabel y mengandung informasi dari x). 

**2. PC merangkum informasi dari actual data sehingga akan dihitung variance-covariance (variance-covariance merepresentasikan informasi yang terkandung pada suatu data)**

```{r}
var(A)
class(var(A))

# alternatif: menggunakan `cov()` untuk menghitung variance-covariance
# cov(A)
# cor(A)
```

Jika diperhatikan variance-covariance yang dihasilkan berupa matriks, dalam aljabar terdapat beberapa operasi perkalian untuk kelas matriks:

- **Perkalian vektor dengan skalar (konstanta)**
```{r}
matrix(c(2,3)) %*% 2
```

Membuat plot untuk melihat vektor awal
```{r}
plot(2, 3)
lines(x = c(0, 2), y = c(0, 3))
```

Membuat plot untuk membandingkan vektor awal dan setelah dilakukan operasi perkalian dengan skalar (konstanta)
```{r}
plot(x = c(2, 4), y = c(3, 6))
lines(x = c(0, 4), y = c(0, 6))
```

> Operasi perkalian antara vektor dengan skalar (konstanta) akan memperbesar vektor dengan arah yang sama.

- **Perkalian matriks dengan vektor**
```{r}
matrix(1:4, nrow = 2) %*% matrix(c(2, 3))
```

Membuat plot untuk membandingkan vektor awal dan setelah dilakukan operasi perkalian dengan matriks
```{r}
plot(x = c(2, 11), y = c(3, 16))
lines(x = c(0, 2), y = c(0, 3))
lines(x = c(0, 11), y = c(0, 16))
```

Jika diperhatikan kedua vektor tersebut seperti berada dalam satu garis karena hampir sejajar, namun jika dicek nilai slope ($y/x$) akan sangat jelas berbeda.
```{r}
# vektor awal
3/2
# vektor setelah dilakukan operasi perkalian dengan matriks
16/11
```

> Operasi perkalian antara matriks dengan vektor akan mengubah besar dan arah vektor

- **Perkalian matriks identitas dengan vektor**
```{r}
matrix(c(1, 0, 0, 1), nrow = 2) %*% matrix(c(2, 3))
```

> Operasi perkalian antara matriks dengan vektor akan mengubah besar dan arah vektor, namun **terdapat matriks yang jika dikalikan dengan vektor tidak akan mengubah besar dan arah vektor, yaitu matriks identitas**.

- **Perkalian matriks rotasi dengan vektor** 
```{r}
matrix(c(-1, 0, 0, -1), nrow = 2) %*% matrix(c(-3, 2))
```

Membuat plot untuk melihat vektor awal
```{r}
plot(-3, 2)
lines(x = c(0, -3), y = c(0, 2))
```

Membuat plot untuk membandingkan vektor awal dan setelah dilakukan operasi perkalian dengan matriks rotasi
```{r}
plot(3, -2)
lines(x = c(0, 3), y = c(0, -2))
```

> Operasi perkalian antara matriks dengan vektor akan mengubah besar dan arah vektor, namun **terdapat matriks yang jika dikalikan dengan vektor hanya akan mengubah arah vektor, yaitu matriks rotasi**.

```{r echo=FALSE, out.width="60%"}
knitr::include_graphics("assets/06-UL/rotation.gif")
```

- **Perkalian matriks dengan eigen vektor**

> Operasi perkalian antara matriks dengan vektor akan mengubah besar dan arah vektor. Namun, terdapat vektor yang unik yang jika dikalikan dengan matriks hasilnya sama dengan mengalikan vektor dengan skalar.

Skalar tersebut adalah konstanta yang disebut sebagai eigen value. Ketika kita mengalikan matriks A dengan eigen vector $x$, hasilnya adalah konstanta $\lambda$ dikalikan $x$. Maka persamaannya adalah:

$$A x = \lambda x$$

- $A$ merupakan matriks variance-covariance 
- $x$ merupakan eigen vector
- $\lambda$ merupkan skalar (konstanta) yang disebut eigen value. 

Eigen vektor tersebut digunakan sebagai transformasi untuk mengubah arah/menentukan arah axis baru yang lebih bisa merangkum informasi pada actual data.

**3. Variance-covariance yang memuat informasi pada actual data akan digunakan untuk memperoleh eigen vektor (transformasi) dan eigen value.**
```{r}
eigen(var(A))
```

Eigen value memuat jumlah informasi yang dirangkum oleh setiap PC, sehingga total eigen value akan sama dengan jumlah variabel pada actual data.
```{r}
1.98271899 + 0.01728101
```

**4. Membuat new data (PC) dari hasil perkalian antara actual data dengan eigen vektor (transformasi).**

New data tersebut akan memuat informasi yang sama dengan actual data dengan arah yang berbeda dari actual data.

$$new\ data = actual\ data \times eigen\ vektor$$

```{r}
# dim(200x2) x dim(2x2)
hasil_pc1 <- A %*% eigen(var(A))$vectors
```