尺度縮減、集群分析、資料視覺化

Exploring Yelp Kaggle Dataset

上週我們對 Yelp Kaggle 裡面的文字做過：

• Bag of Words
• Sentiment Analysis
• Linguistic Inquiry & Word Count (by Empath Text Classifier, LIWC alike)

這一些 文字分析 之後，本周我們繼續用這一組資料，來示範：

• 尺度縮減
• 集群分析
• 資料視覺化

這幾種分析技術的綜合運用。

---

library(magrittr)
library(Rtsne)
library(RColorBrewer)
library(randomcoloR)
library(wordcloud)
library(d3heatmap)
library(igraph)  
library(reshape2)
library(highcharter)

load('data/yelp1.rdata')  # loading yelp data & sentiment scores
load('data/empath.rdata') # loading empath scores

載入 packages 和 data 之後，一開始我們有 6 個 Data Frame:

Data Frame	No. Rows x Cols	Objects
biz	11,537 x 11	businesses
user	43,873 x 7	users
check	262,764 x 4	check-in’s
review	215,879 x 8	reviews
senti	215,879 x 10	reviews’ sentiment scores
scores	215,879 x 194	reviews’ Empath/LIWC scores

1. 商店的類別 (508 Business Categories)

這份資料裡面有11,537家商店(Business)和508種商業類別(Categories)， 一家商店可以同時隸屬於很多(0 ~ 10)個類別， 我們首先考慮：

• 每一個商業類別裡面有多少家商店？ 這些商店總共被評論過多少？
• 這些商業類別之間有相關性嗎？ 哪一些商業類別經常會同時出現？

(1a) 商業類別資料整理

# number of biz per category
CA = biz$cat %>% strsplit('|',T) %>% unlist %>% table %>% 
  data.frame %>% 'names<-'(c('name','nbiz'))
# number of review per category
CA$nrev = CA$name %>% sapply(function(z){sum(
  review$bid %in% biz$bid[grep(z,biz$cat,fixed=T)] )})
# average number of reviews per business
CA$avg.rev = CA$nrev / CA$nbiz    
CA = CA[order(-CA$nrev),]  # order CA by no. review
rownames(CA)= CA$name 
CA$name = NULL

從biz$categoey裡面整理出 508 個 Categories， 並算出每一個 Category 的：

• nbiz : number of business (in the category)
• nrev : number of revierw (in the category)
• avg.rev : average number of reviews per business

放在 CA 這個 Data Frame 裡面：

CA

# category-business matrix
mxBC = rownames(CA) %>% sapply(function(z)
  grepl(z,biz$cat,fixed=T))
rownames(mxBC) = biz$bid
dim(mxBC) 

[1] 11537   508

將 Business(11,537) 和 Category(508) 的對應關係放在mxBC裡面。

(1b) 尺度縮減 (Dimension Reduction)

使用tSNE，將mxBC的尺度 [11,537 x 508] 縮減為 [2 x 508] …

t0 = Sys.time()
set.seed(123)
tsneCat = ifelse(mxBC,1,0) %>% t %>% 
  Rtsne(check_duplicates=F,theta=0.0,max_iter=3000)
Sys.time() - t0

Time difference of 23.86 secs

(1c) 階層式集群分析 (Hierarchical Clustering)

在縮減尺度之中做階層式集群分析。

Y = tsneCat$Y           # tSNE coordinates
d = dist(Y)             # distance matrix
hc = hclust(d)          # hi-clustering

K = 60                  # number of clusters 
CA$grp = cutree(hc,K)   
table(CA$grp)

 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 
 5  2 10  7  9  6  4  6 14  8  3  2  5  8  9  3  2  7 11 12 10  3 15  5  4  6  7  8  7 
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 
 9  7 12  6  5 11  9 11 21 18 12  8 10  8 14 15  2  5 11  6 12 15  6  4  8  4 13 15 13 
59 60 
10 10 

經重覆嘗試之後，將這508個商業類別(Categories) 分成60個商業類別群組(Category Groups)。

(1d) 字雲 (Word Cloud)

pals = distinctColorPalette(K)  # palette of K colors 
png("category.png", width=3200, height=1800)
textplot(Y[,1],Y[,2],rownames(CA),font=2, 
         col = pals[CA$grp],         # color by group    
         cex = (0.5+log(CA$nrev)/5)) # size by no. reviews
dev.off()

null device 
          1 

將字雲畫在category.png裡面：

• 每個字代表一個商業類別(Categories)
• 字的顏色代表商業類別群組(Category Groups)
• 字的大小代表這個商業類別被評論的次數 (number of reviews)
• 靠在一起的、同一種顏色的字，代表經常一起出現的商業類別

.

2. 評論的內容 (194 Content Classes)

接下來考慮評論的內容，上週我們已經使用 Empath Text Classifier ，依其預設的194種內容(Class)， 對這215,879篇評論分別做過評分。 也就是說，文集之中的每一篇評論都有194個內容評分， 存放在scores這個Data Frame裡面。

dim(scores)

[1] 215879    194

使用這一些資料，我們可以對這194種內容(Classes)進行分群， 也可以用字雲來呈現不同內容之間的相關性。

(2a) 內容權重 (Class Weights)

計算每一種內容的權重(Weight: the sum of class scores within the corpus)， 放在wClass裡面，並將scores(內容評分資料)依權重排序。

# 
# order the score matrix by class weights
scores = scores[,order(-colSums(scores))]
wClass = colSums(scores)  # class weights
head(wClass,20)

          eating          cooking       restaurant         shopping positive_emotion 
          7289.6           6286.9           6121.8           3248.1           2733.0 
         friends         business           giving            party         vacation 
          2693.6           1996.2           1780.7           1726.1           1642.5 
        optimism      achievement   shape_and_size negative_emotion       occupation 
          1548.4           1457.1           1382.8           1370.4           1350.9 
     celebration        traveling             home         children           family 
          1248.4           1235.9           1127.3           1106.8           1098.7 

(2b) 尺度縮減 (Dimension Reduction)

使用tSNE，將scores的尺度 [216879 x 194] 縮減為 [2 x 194] …

t0 = Sys.time()
set.seed(123)
tsneClass = scores %>% scale %>% as.matrix %>% t %>% 
  Rtsne(theta=0.0,max_iter=3000)
Sys.time() - t0

Time difference of 46.549 secs

(2c) 階層式集群分析 (Hierarchical Clustering)

在縮減尺度之中做階層式集群分析。

Y = tsneClass$Y      # tSNE coordinates
d = dist(Y)          # distance matrix
hc = hclust(d)       # hi-clustering

K = 40
gpClass = cutree(hc,K)  # K groups
table(gpClass)          # no. classes per group

gpClass
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 
 3  9  6  5  4  7  4  3  8  6  4  3  5  6 10  4  4  3  4  6  8  9  2  3  7  5  4  6  4 
30 31 32 33 34 35 36 37 38 39 40 
 4  3  2  4  7  2  5  4  3  5  3 

返覆嘗試之後，將這194種內容(Classes) 分成40個內容群組(Class Groups)。

(2d) 字雲 (Word Cloud)

pals = distinctColorPalette(K)  # palette of K colors 
png("classes.png", width=3200, height=1800)
textplot(Y[,1],Y[,2],colnames(scores),font=2, 
  col = pals[gpClass],      # color by group    
  cex = 0.5+log(wClass)/3)  # size by class weight
dev.off()

null device 
          1 

將字雲畫在classes.png裡面：

• 每個字代表一種內容(Class)
• 字的顏色代表內容群組(Class Groups)
• 字的大小代表這種內容在文集之中的權重
• 靠在一起的、同一種顏色的字，代表經常一起出現的內容

.

(2e) 內容之間的相關性

繼續分析之前，先用熱圖檢查一下內容評分之間的相關性 …

# correlation between the classes
cr = cor(scores)  
cr %>% d3heatmap(scale='none',col=cm.colors(256))

找出相關係數較高的內容種類

corgrp = function(x, threshold=0.6) {
  x = x*lower.tri(x)
  check = which(x>threshold,arr.ind=T)
  gcr = graph.data.frame(check, directed=F)
  gcr = split(unique(as.vector(check)),clusters(gcr)$membership)
  lapply(gcr, function(g){rownames(x)[g]})  }
corgrp(cr, 0.6) 

$`1`
[1] "cooking"    "restaurant" "eating"    

$`2`
[1] "achievement"      "positive_emotion"

$`3`
[1] "giving"     "occupation" "phone"      "business"  

$`4`
[1] "celebration" "party"      

$`5`
[1] "affection" "optimism"  "love"     

$`6`
[1] "payment"   "valuable"  "economics" "money"    

$`7`
[1] "hygiene"  "cleaning"

$`8`
[1] "vehicle" "car"     "driving"

$`9`
[1] "pain"           "shame"          "suffering"      "swearing_terms"
[5] "violence"       "emotional"      "hate"          

$`10`
[1] "hearing" "noise"   "sound"   "listen" 

$`11`
[1] "toy"  "play"

$`12`
[1] "leader" "order" 

$`13`
[1] "dance"   "musical" "music"  

$`14`
[1] "fire"   "warmth"

$`15`
[1] "water"    "swimming" "sailing"  "exotic"   "ocean"   

$`16`
[1] "technology"  "programming" "computer"    "internet"   

$`17`
[1] "school"  "college"

$`18`
[1] "nervousness" "contentment"

$`19`
[1] "disgust" "anger"  

3. 商業類別 與 評論內容

接下來我們可以做商業類別(508 Categories)和評論內容(194 Classes) 之間的交叉分析。

(3a) 交叉查詢範例

先用以下這一個範例來示範交叉查詢的做法

• 在評論數超過500的商業類別之中，哪些類別的評論之中的正面情緒和負面情緒的相關係數是最高的呢？

# correlation between positive & negative emotion 
mxBC[,CA$nrev>500] %>% apply(2,function(v) {
  i = review$bid %in% rownames(mxBC)[v]
  cor(scores$positive_emotion[i], scores$negative_emotion[i])
}) %>% sort %>% tail(20)

                Pakistani            Middle Eastern                 Soul Food 
                 0.035701                  0.041726                  0.042036 
               Pet Stores                  Day Spas                  Car Wash 
                 0.043151                  0.051408                  0.053279 
               Drugstores                   Massage                Automotive 
                 0.055246                  0.078279                  0.089250 
                    Cafes               Auto Repair         Stadiums & Arenas 
                 0.091302                  0.097786                  0.102325 
Cosmetics & Beauty Supply                 Skin Care                  Dentists 
                 0.105822                  0.108231                  0.138033 
                    Tires          Health & Medical              Pet Services 
                 0.145028                  0.182872                  0.218834 
                     Pets                   Doctors 
                 0.244424                  0.275944 

我們可以看到Doctors和Pets這兩種商業類別的評論最常會同時出現正面和負面情緒。

(3b) 各商業類別的內容權重

不同商業類別的評論裡面會有不一樣的內容， 我們可以用熱圖來呈現各商業類別的內容分布狀況。 先把內容評分依商業類別平均起來， 放在wxClass這個矩陣裡面。

wxClass = apply(mxBC,2,function(v){     # for every category
  i = review$bid %in% rownames(mxBC)[v] # find its reviews 
  colMeans(scores[i,])                  # average their class scores
})
dim(wxClass)

[1] 194 508

由於wxClass這個矩陣太大， 我們只畫出評論數最多的50個商業類別和權重最大的50種評論內容

wxClass[1:50,1:50] %>% log %>% 
  d3heatmap(T,T,scale='none',col='PiYG')

(3c) 群組熱圖

因為商業類別(508 categories)和評論內容(194 classes)的數量都很多 (這就是大數據的特徵)， 我們要在群組這個層面，才比較容易觀察整個文集的內容分布狀況。 其實，這就是大數據分析之中，我們常常需要先做集群分析的理由。 我們用以下的熱圖呈現商業類別群組(60)和內容群組(40)之間的關係。

x = matrix(0, nrow=max(gpClass), ncol=max(CA$grp),
  dimnames=list(sprintf('CLS%02d',1:max(gpClass)),
                sprintf('CAT%02d',1:max(CA$grp))))
for(i in 1:nrow(x)) for(j in 1:ncol(x))  
  x[i,j] = sum( wxClass[gpClass==i, CA$grp==j] )
t(x) %>% {log(0.005+.)} %>% d3heatmap(scale='none',col='PiYG')

從以上的群組熱圖裡面我們可以清楚的看到整個文集(在各商業類別之中)的內容分布狀況。
如果需要看某一群組之中有哪一些商業類別或評論內容，可以這樣做：

rownames(CA)[CA$grp==8]

[1] "Sandwiches"             "Delis"                  "Bagels"                
[4] "Sporting Goods"         "Bikes"                  "Food Delivery Services"

names(scores)[gpClass==1]

[1] "eating"     "cooking"    "restaurant"

4. 趨勢分析

(4a) Trend of Content Classes

txClass = tx = split(scores, cut(review$date,'quarter')) %>% 
  sapply(colSums) %>% apply(2, function(v) v/sum(v))
# Trend of Classes
df = data.frame(class=rownames(tx), tx[,9:32]) %>% melt('class')
df$date = as.Date(
  substr(as.character(df$variable),2,11),format="%Y.%m.%d")
df$value = round(100*df$value,3)
subset(df, class %in% rownames(tx)[1:194]) %>% 
  hchart("spline",hcaes(x=date,y=value,group=class)) %>% 
  hc_legend(align="left", layout="vertical",verticalAlign="top") %>% 
  hc_add_theme(hc_theme_flat()) %>% 
  hc_title(text='Weights of Classes (% of Total Weight)')

(4b) Trend of Class Groups

# Trend of Class Group
x = split(data.frame(tx),gpClass) %>% sapply(colSums) %>% t
df = data.frame(
  class=sprintf('G%02d',1:max(gpClass)), x[,9:32]) %>% 
  melt('class')
df$date = as.Date(
  substr(as.character(df$variable),2,11),format="%Y.%m.%d")
df$value = round(100*df$value,3)
df %>% hchart("spline",hcaes(x=date,y=value,group=class)) %>% 
  hc_legend(align="left", layout="vertical",verticalAlign="top") %>% 
  hc_add_theme(hc_theme_flat()) %>% 
  hc_title(text='Weights of Class Groups (% of Total Weight)')

(4c) Trend of Business Categories

txCat = tx = split(review, cut(review$date,'quarter')) %>% 
  sapply(function(x) apply(mxBC,2,function(v)
      sum(x$bid %in% rownames(mxBC)[v]) )) %>% 
  apply(2, function(v) v/sum(v))
# Trend of Categories
df = data.frame(category=rownames(tx), tx[,9:32]) %>% melt('category')
df$date = as.Date(
  substr(as.character(df$variable),2,11),format="%Y.%m.%d")
df$value = round(100*df$value,3)
subset(df, category %in% rownames(CA)[1:30]) %>% 
  hchart("spline",hcaes(x=date,y=value,group=category)) %>% 
  hc_legend(align="left", layout="vertical",verticalAlign="top") %>% 
  hc_add_theme(hc_theme_flat()) %>% 
  hc_title(text='Weights of Categories (% of Total Reviews)')

(4d) Trend of Business Category Groups

# Trend of Category Groups
x = split(data.frame(tx),CA$grp) %>% sapply(colSums) %>% t
df = data.frame(
  category=sprintf('G%02d',1:max(CA$grp)), x[,9:32]) %>% 
  melt('category')
df$date = as.Date(
  substr(as.character(df$variable),2,11),format="%Y.%m.%d")
df$value = round(100*df$value,3)
df %>% hchart("spline",hcaes(x=date,y=value,group=category)) %>% 
  hc_legend(align="left", layout="vertical",verticalAlign="top") %>% 
  hc_add_theme(hc_theme_flat()) %>% 
  hc_title(text='Weights of Category Groups (% of Total Reviews)')