AS6-0 集群分析

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

[1] "C"

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('data/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
plot(hc)
rect.hclust(hc, k=5, border="red")

W$group = cutree(hc, k=5) %>% factor
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
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("data/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)

              Title Unknown Action Adventure Animation Childrens Comedy Crime Documentary
1  Toy Story (1995)       0      0         0         1         1      1     0           0
2  GoldenEye (1995)       0      1         1         0         0      0     0           0
3 Four Rooms (1995)       0      0         0         0         0      0     0           0
4 Get Shorty (1995)       0      1         0         0         0      1     0           0
5    Copycat (1995)       0      0         0         0         0      0     1           0
  Drama Fantasy FilmNoir Horror Musical Mystery Romance SciFi Thriller War Western
1     0       0        0      0       0       0       0     0        0   0       0
2     0       0        0      0       0       0       0     0        1   0       0
3     0       0        0      0       0       0       0     0        1   0       0
4     1       0        0      0       0       0       0     0        0   0       0
5     1       0        0      0       0       0       0     0        1   0       0

sum(M$Comedy)             # 喜劇片

[1] 502

sum(M$Western)            # 西部片

[1] 27

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

[1] 863

1.3 距離矩陣

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)
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 

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

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)

【問題討論】

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

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

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

---

2. Flower Image

2.1 整理資料

# Read data
flower = read.csv("data/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("data/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: cannot allocate vector of size 498.0 Gb

【Q】 What is the problem?

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("data/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)        # 建立模型
tumorClusters = predict(KMC.kcca, newdata = tumorVector) # 進行預測(轉換)
Sys.time() - t0

Time difference of 31.98 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集群分析有什麼差異？ 它們分別用在什麼狀況？

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

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

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

---