AS6-0 集群分析

Sys.setlocale("LC_ALL","C")

[1] "C/C/C/C/C/en_US.UTF-8"

packages = c(
  "dplyr","ggplot2","d3heatmap","googleVis","devtools","plotly", "xgboost",
  "magrittr","caTools","ROCR","corrplot", "rpart", "rpart.plot",
  "doParallel", "caret", "glmnet", "Matrix", "e1071", "randomForest",
  "flexclust", "FactoMineR", "factoextra"
  )
existing = as.character(installed.packages()[,1])
for(pkg in packages[!(packages %in% existing)]) install.packages(pkg)

rm(list=ls(all=T))
options(digits=4, scipen=12)
library(dplyr)
library(ggplot2)
library(flexclust)
library(FactoMineR)
library(factoextra)

---

A. 集群分析與尺度縮減

A1. 批發交易資料

W = read.csv("wholesales.csv")
W$Channel = factor( paste0("Ch",W$Channel) )
W$Region = factor( paste0("Reg",W$Region) )
W[3:8] = lapply(W[3:6], log, base=10)
summary(W)

 Channel    Region        Fresh            Milk         Grocery          Frozen    
 Ch1:298   Reg1: 77   Min.   :0.477   Min.   :1.74   Min.   :0.477   Min.   :1.40  
 Ch2:142   Reg2: 47   1st Qu.:3.495   1st Qu.:3.19   1st Qu.:3.333   1st Qu.:2.87  
           Reg3:316   Median :3.930   Median :3.56   Median :3.677   Median :3.18  
                      Mean   :3.792   Mean   :3.53   Mean   :3.666   Mean   :3.17  
                      3rd Qu.:4.229   3rd Qu.:3.86   3rd Qu.:4.028   3rd Qu.:3.55  
                      Max.   :5.050   Max.   :4.87   Max.   :4.967   Max.   :4.78  
 Detergents_Paper   Delicassen  
 Min.   :0.477    Min.   :1.74  
 1st Qu.:3.495    1st Qu.:3.19  
 Median :3.930    Median :3.56  
 Mean   :3.792    Mean   :3.53  
 3rd Qu.:4.229    3rd Qu.:3.86  
 Max.   :5.050    Max.   :4.87  

A2. 兩個區隔變數

hc = W[,3:4] %>% scale %>% dist %>% hclust  
# [,3:4]: fresh, milk
# scale: Normalize 標準化（mean=0, sd=1)
# dist: 距離矩陣
# hclust: 集群分析（h）
# 回傳一個 Hierarchical Clustering 的物件
plot(hc)
rect.hclust(hc, k=5, border="red")

# k 是測試出來的，可以多方測試看切得好不好再來決定～
# 找切下去「垂直距離較長的」通常比較好～
# 沒有標準知道切多少比較好

W$group = cutree(hc, k=5) %>% factor
# cutree: 分群
ggplot(W, aes(x=Fresh, y=Milk, col=group)) +
  geom_point(size=3, alpha=0.5) + 
  theme_light()

A3. 六個區隔變數

hc = W[,3:7] %>% scale %>% dist %>% hclust
# [,3:7]: Fresh,......, Detergent_Paper
plot(hc)
W$group = factor(cutree(hc, k=8))
rect.hclust(hc, k=8, border="red")

library(FactoMineR)
library(factoextra)
fviz_dend(
  hc, k=8, show_labels=F, rect=T, rect_fill=T,
  labels_track_height=0,
  palette="ucscgb", rect_border="ucscgb")

A4. 尺度縮減

Dimension Reduction with PCA (Principle Component Analysis, 主成分分析)

W[,3:8] %>% PCA(graph=F) %>% fviz_pca_biplot(
  label="var", col.ind=W$group,
  pointshape=19, mean.point=F,
  addEllipses=T, ellipse.level=0.7,
  ellipse.type = "convex", palette="ucscgb",
  repel=T
  )

---

1. Cluster Analysis for Movies

主要議題：依類型(Genre)對電影分類

學習重點：

• 集群分析的基本觀念
• 距離矩陣：Distance Matrix
• 層級式集群分析：Hierarchical Cluster Analysis
• 樹狀圖(Dendrogram)的判讀
• 依據樹狀圖決定要分多少群
• 以群組平均值檢視各族群的屬性

1.1 整理資料

M = read.table("movieLens.txt", header=FALSE, sep="|",quote="\"")
# Assign column names
colnames(M) = c(
  "ID", "Title", "ReleaseDate", "VideoReleaseDate", "IMDB", 
  "Unknown", "Action", "Adventure", "Animation", "Childrens", 
  "Comedy", "Crime", "Documentary", "Drama", "Fantasy", "FilmNoir", 
  "Horror", "Musical", "Mystery", "Romance", "SciFi", "Thriller",
  "War", "Western")
# Remove unnecessary variables
M$ID = NULL
M$ReleaseDate = NULL
M$VideoReleaseDate = NULL
M$IMDB = NULL
# Remove duplicates
M = unique(M)

1.2 檢視資料

head(M, 5)

sum(M$Comedy)             # 喜劇片

[1] 502

sum(M$Western)            # 西部片

[1] 27

sum(M$Romance | M$Drama)  # 浪漫劇情片

[1] 863

# ｜： or

1.3 距離矩陣

# 沒有scale: 因為結果都是0、1，沒有必要做 normalize
dmx= dist(M[2:20], method="euclidean")
dmx %>% as.matrix %>% dim

[1] 1664 1664

1.4 層級式集群分析

hclust1 = hclust(dmx, method = "ward.D") 

1.5 檢視樹狀圖

plot(hclust1)
rect.hclust(hclust1, k=5, border="red")

1.6 切割群組

grp = cutree(hclust1, k = 5)
# cutree 出來的結果會是「整數向量」
table(grp)

grp
  1   2   3   4   5 
824 370 209 196  65 

1.7 檢查群組屬性

# 做完「集群分析」後會想做的第一件事
tapply(M$Action, grp, mean)

      1       2       3       4       5 
0.28641 0.00000 0.00000 0.06633 0.00000 

# 統計族群中【Action群】的比例
tapply(M$Romance, grp, mean)

      1       2       3       4       5 
0.05825 0.00000 0.00000 1.00000 0.00000 

1.8 The sapply-split-... Combo：

sapply(split(M[,2:20], grp), colMeans) %>% round(3)

                1 2 3     4 5
Unknown     0.002 0 0 0.000 0
Action      0.286 0 0 0.066 0
Adventure   0.161 0 0 0.000 0
Animation   0.051 0 0 0.000 0
Childrens   0.146 0 0 0.000 0
Comedy      0.177 0 1 0.418 1
Crime       0.123 0 0 0.031 0
Documentary 0.061 0 0 0.000 0
Drama       0.238 1 0 0.434 1
Fantasy     0.027 0 0 0.000 0
FilmNoir    0.028 0 0 0.005 0
Horror      0.107 0 0 0.010 0
Musical     0.068 0 0 0.000 0
Mystery     0.073 0 0 0.000 0
Romance     0.058 0 0 1.000 0
SciFi       0.121 0 0 0.000 0
Thriller    0.279 0 0 0.092 0
War         0.086 0 0 0.000 0
Western     0.033 0 0 0.000 0

# [,2:20]: 去掉片名的19個標籤
# sapply（集合, ）: 把第一個list中的每個項目都去呼叫後面的 grp
# (split(M[]), grp): 把一個資料框根據 grp 切開，分成5個（k分的個數）
# colMean: 同 tapply(M$Action, grp, mean)，共做了 5 次
# round(3): 小數點後三位

1.9 資料視覺化

layout(matrix(c(1,2,2), 3, 1))
par(mar=c(2,3,1,1), cex=0.8)
table(grp) %>% barplot(col=3:7, names.arg=paste0("Group-",1:5))
par(mar=c(6,3,2,1))

sapply(split(M[,2:20], grp), colMeans) %>% t %>% 
  barplot(beside=T, col=3:7, las=2)

【問題討論】

從管理的角度來看，我們為甚麼要分群？

• 管理面：藉由分群的方式，讓相同類型或屬性的人員放在一起，以便於管理者進行管理
• 成本面：當把相類似的群體放咋一起時，單一管理者可以同時管理最多人，以縮減人事成本
• 執行面：藉由先將類似預測對象分群針對某一群人進行分群，結果通常會變得更加準確

我們為甚麼要做尺度縮減？

• 將資料投射在變數平面，讓人們得以辨識「不同群體」，再將n-dimension的資料壓縮到2D平面，使得會壓在一個資料損失最少的角度，讓2D資料呈現量最大。

我們要如何把集群分析的結果轉化為策略呢？

• 藉由集群分析進行尺度縮減之後，可以將n-D資料壓縮2在資料呈現量最大的2D平面上。
• 壓縮後的2D圖型，可以可以更加視覺化的讓我們瞭解資料彼此之間的關聯性。
• 藉由瞭解資料彼此間的關聯性，可以藉由資料彼此間的關係，或是資料與時間的關聯性，來針對使用者投放個人化的行銷策略，以增加銷售額、毛利、品牌印象，得以讓策略

---

2. Flower Image

2.1 整理資料

# Read data
flower = read.csv("flower.csv", header=FALSE)
# Change the data type to matrix
flowerMatrix = as.matrix(flower)
dim(flowerMatrix)

[1] 50 50

# Turn matrix into a vector
flowerVector = as.vector(flowerMatrix)
length(flowerVector)

[1] 2500

2.2 距離矩陣

# Compute distances
distance = dist(flowerVector, method = "euclidean")

2.3 層級式集群分析

# Hierarchical clustering
clusterIntensity = hclust(distance, method="ward.D")

2.4 樹狀圖

# Plot the dendrogram
plot(clusterIntensity)
# Select 3 clusters
rect.hclust(clusterIntensity, k = 3, border = "red")

切割群組

flowerClusters = cutree(clusterIntensity, k = 3)
table(flowerClusters)

flowerClusters
   1    2    3 
1634  272  594 

# flowerClusters

族群平均(畫素顏色深淺度)

# Find mean intensity values
tapply(flowerVector, flowerClusters, mean)

      1       2       3 
0.08574 0.50826 0.93148 

圖像比較

# Plot the image and the clusters
dim(flowerClusters) = c(50,50)
par(mfrow=c(1,2), mar=c(2,2,2,2))
# Original image
image(flowerMatrix,axes=FALSE,col=grey(seq(0,1,length=256)),main="Original")
# New image
image(flowerClusters, axes = FALSE, main="3 Cluster")

---

3. MRI Image

3.1 整理資料

# Read data
healthy = read.csv("healthy.csv", header=FALSE)
healthyMatrix = as.matrix(healthy)
dim(healthyMatrix)

[1] 566 646

3.2 畫出圖形

# Plot image
par(mar=c(1,1,1,1))
image(healthyMatrix,axes=FALSE,col=grey(seq(0,1,length=256)))

3.3 距離矩陣

# Compute distances
healthyVector = as.vector(healthyMatrix)
distance = dist(healthyVector, method = "euclidean")

Error: vector memory exhausted (limit reached?)

【Q】 What is the problem?

• 圖形的pixel多，矩陣乘出來的數值會變得非常非常大，因而導致出現error的訊息以防止電腦當機。

3.4 KMeans集群分析

# Run k-means
k = 5
set.seed(1)
KMC = kmeans(healthyVector, centers = k, iter.max = 1000)

3.5 檢查分群結果

# View(KMC)
table(KMC$cluster)

     1      2      3      4      5 
 20556 101085 133162  31555  79278 

KMC$centers

     [,1]
1 0.48177
2 0.10619
3 0.01962
4 0.30943
5 0.18421

3.6 畫出分群結果

# Extract clusters
X = KMC$cluster
# Plot the image with the clusters
dim(X) = c(nrow(healthyMatrix), ncol(healthyMatrix))
# Plot image
par(mar=c(1,1,1,1))
image(X, axes = FALSE, col=rainbow(k))

3.7 讀進、轉換測試圖形

tumor = read.csv("tumor.csv", header=FALSE)
tumorMatrix = as.matrix(tumor)
dim(tumorMatrix)

[1] 571 512

tumorVector = as.vector(tumorMatrix)
length(tumorVector)

[1] 292352

3.8 將原圖形之分群規則套用到測試圖形

# Apply clusters from before to new image, using the flexclust package
library(flexclust)
t0 = Sys.time()
KMC.kcca = flexclust::as.kcca(KMC, healthyVector)        # 建立模型

Found more than one class "kcca" in cache; using the first, from namespace 'flexclust'
Also defined by 'kernlab'
Found more than one class "kcca" in cache; using the first, from namespace 'flexclust'
Also defined by 'kernlab'

tumorClusters = predict(KMC.kcca, newdata = tumorVector) # 進行預測(轉換)

Found more than one class "kcca" in cache; using the first, from namespace 'flexclust'
Also defined by 'kernlab'

Sys.time() - t0

Time difference of 34.51 secs

3.9 圖像比較

# Visualize the clusters
dim(tumorClusters) = c(nrow(tumorMatrix), ncol(tumorMatrix))
par(mfrow=c(1,2), mar=c(1,1,2,1))
image(X, axes = FALSE, col=rainbow(k), main="Healthy")
image(t(tumorClusters)[,571:1], axes = FALSE, col=rainbow(k), main="Tumor")

【學習重點】

• 集群分析在圖像處理的應用
• 單區隔變數的集群分析
• 集群分析模型

【問題討論】

層級式和K-Means集群分析有什麼差異？ 它們分別用在什麼狀況？

• 層級式：
  • 從單筆資料開始，不用一開始就決定要分成多少群，最終所有的clusters都會成為一個cluster（通常適用於不要太大的資料集）
• K-Means：
  • 事先定義好切成k組，最終即會出現所需分組組數的分群結果即便（資料量再大都能輕易做出來 ）

集群分析模型和普通的集群分析有什麼差異？

• 集群分析模型：

• 集群分析模型與普通集群分析模型不同，集群分析模型可以將原始資料也分成Training Data與Testing Data來建立模型，再由建立出來的集群分析模型將更多未知的資料為進去，進而能夠針對不同群族進行結果之預測

• 普通集群分析：
  • 目標由資料分群結果而定，而非來自於training data預測結果之模型，藉由兩類型（層級式、K-means）的建模方式，可以將資料依照階層將所有資料進行分層繪製出dendorgram，或是依照欲分組數k得出分類結果

什麼時候需要建集群分析模型？ 集群分析模型的用法？

• 當資料當中族群組內差異甚大，但是這其中又有一群一群的關係，使分組後組內差異小的情況產生時，便適合使用建置集群分析模型，若是與其他監督式學習方法一同執行，所做出預測之準確度會遠比沒有建集群分析模型高出許多，進而也能幫助做出決策分析。

圖像處理和圖像辨識有什麼差異？

• 圖像處理：
  • 圖像處理是較基礎的圖像製作，藉由將圖像們細分成不同的pixel後組合而成，可以讓所有的圖像都能被分割成不同的小單位，且組成的每一單位佔原圖像的比例也攸關到圖像的解析度。（每一張圖像的內部處理）
• 圖像辨識：
  • 圖像辨識是機器學習所採用的概念，藉由將許許多多類似的圖片丟進辨識的工具後，能夠辨別出形狀、顏色等細節的差異，而若是丟足夠多的圖像進去機器的話，其辨識準確度會非常高。（多張圖像的疊合對比處理）

---