👨‍🏫

點我看教學影片 (ctrl + click)


rm(list=ls(all=T))
options(digits=4, scipen=40)
library(dplyr)
Load Yelp10 Data
load("data/Biz.rdata")
load("data/Rev.rdata")
LOAD = TRUE
if(LOAD) { load("data/tsne.rdata") }

(1) rev - 評論資料框

評論資料框的表頭 + 評論資料框包括評論id、商店id、評論星等、評論時間、評論情緒等等資訊。

head(rev)
  rid    bid       date stars cool funny useful anger anticipation disgust fear joy sadness
1   1 154938 2014-02-17     4    0     0      1     1            1       0    0   3       0
2   2  36983 2016-07-24     4    0     0      0     0            1       1    0   1       0
3   3 103253 2012-04-07     5    0     0      2     0            1       0    0   1       0
4   4  25101 2015-09-11     5    0     0      0     0            1       1    0   0       0
5   5 182709 2016-01-22     5    1     1      1     0            1       0    0   1       0
6   6 135863 2014-09-17     5    0     0      1     0            2       1    0   6       0
  surprise trust negative positive            business_id                user_id nchar year
1        0     3        1        3 Ue6-WhXvI-_1xUIuapl0zQ gVmUR8rqUFdbSeZbsg6z_w   322 2014
2        1     1        1        3 Ae4ABFarGMaI5lk1i98A0w Y6qylbHq8QJmaCRSlKdIog   137 2017
3        1     1        0        1 lKq4Qsz13FDcAVgp49uukQ SnXZkRN9Yf060pNTk1HMDg   108 2012
4        0     0        1        1 6nKR80xEGHYf2UxAe_Cu_g VcmSgvslHAhqWoEn16wjjw   113 2016
5        0     1        0        2 Z_mJYg3vi8cPZHa1J4BALw NKF9v-r0jd1p0JVi9h2T1w   200 2016
6        1     5        0        7 R1PQEK6qvrZVC9qcWfKvDA c2MQ_LPuvtiiKFR_-OY9pg   326 2015
1.1 Quick Check
  • 評論人數:1518169個使用者
n_distinct(rev$user_id)      # no. user = 1518169
[1] 1518169
  • 商店數量:188593間商店
n_distinct(rev$bid)          # no. biz =   188593
[1] 188593
  • 查看每則評論的長度
sapply(c(5, 10, 50, 100), function(i) sum(rev$nchar <= i))
[1]     69    163   8630 177993
breaks = as.Date(c("2004-07-02", paste0(2005:2018, "-07-02"))) # 以07/02這一天來切分年份
rev$year = as.integer(cut(rev$date, breaks)) + 2004 # 新增year欄位
  • 以年份劃分評論總數
  • 以各個星等劃分評論總數
  • 評論回應字數的分布
par(cex=0.8, mfrow=c(1,3), mar=c(7,5,4,2))
table(rev$year) %>% 
  barplot(las=2, main="#Reviews by Year(cut at Jul02)", 
          xlab="", ylab="")
table(rev$stars) %>% barplot(main="No. Stars")
hist(rev$nchar, main="No.Characters")

  • 以下分別是將stars星等、cool酷、funny有趣、useful實用取其平均分數
# average scores
df = aggregate(cbind(stars,cool,funny,useful) ~  year, data = rev, FUN = mean)
par(cex=0.8, mfrow=c(1,4), mar=c(3,4,4,1))
mapply(barplot, df[2:5], main=names(df)[2:5], las=2)

      stars cool funny useful
 [1,]   0.7  0.7   0.7    0.7
 [2,]   1.9  1.9   1.9    1.9
 [3,]   3.1  3.1   3.1    3.1
 [4,]   4.3  4.3   4.3    4.3
 [5,]   5.5  5.5   5.5    5.5
 [6,]   6.7  6.7   6.7    6.7
 [7,]   7.9  7.9   7.9    7.9
 [8,]   9.1  9.1   9.1    9.1
 [9,]  10.3 10.3  10.3   10.3
[10,]  11.5 11.5  11.5   11.5
[11,]  12.7 12.7  12.7   12.7
[12,]  13.9 13.9  13.9   13.9
[13,]  15.1 15.1  15.1   15.1
[14,]  16.3 16.3  16.3   16.3

(2) user - 評論人資料框

  • 根據每個user_id，計算其所有回覆的評論之評論數量、平均星等、平均funny的分數、平均useful的分數、平均cool的分數
user = rev %>% group_by(user_id) %>% summarise(
  n = n(),
  star = mean(stars),
  funny = mean(funny),
  useful = mean(useful),
  cool = mean(useful)
  )
save(user, file="data/User.rdata", compress=T)
  • 以圖示來呈現上述使用者評論相關之各個頻率分布
par(cex=0.8, mfrow=c(1,2), mar=c(5,5,4,2))
hist(log(user$n), main="No. Reviews per User (log)") # y軸為頻率
hist(user$star, main="Avg. Stars per User (log)")

par(cex=0.8, mfrow=c(1,3), mar=c(5,5,4,2))
hist(pmin(user$funny,10), main="Avg. Funny's per User")
hist(pmin(user$cool,10), main="Avg. Cool's per User")
hist(pmin(user$useful,10), main="Avg. Useful's per User")

(3) X - 商店類別矩陣 Biz-Category Matrix

3.1 X - BC matrix, 1306 categories
  • 每一個商店可能屬於很多個商業類別，所以商店和類別之間的關係需要用矩證的方式表示。
dim(X)
Loading required package: Matrix
NULL
  • 如下圖所示，大多數的商店都屬於多個類別
par(cex=0.8, mar=c(3,4,4,2))
rowSums(X) %>% table %>% head(10) %>% barplot(main="No. Categoy per Biz")

  • 各類別的商店數大致上是長尾分佈(power distribution)
par0 = par(cex=0.7, mar=c(11,4.5,3,0))
colSums(X)[1:40] %>% barplot(las=2, main="Top 40 Biz Category")

3.2 X - dense BC matrix, 936 categories
  • 有一些商業類別的商店很少，我們決定只留下商店數大於20的商業類別
  • 188593間商店，936種商業類別
X = X[,colSums(X) > 20]
dim(X)                   # 188593    936
[1] 188593    936
identical(B$business_id, rownames(X))  # TRUE
[1] TRUE
3.3 C - 商業類別摘要
C = apply(X, 2, function(v) c(sum(v), sum(B[v,]$review_count)))
C = C %>% t %>% data.frame %>% setNames(c("n_biz", "n_rev")) %>% 
  mutate(a_rev = n_rev/n_biz)
Warning: `as_dictionary()` is soft-deprecated as of rlang 0.3.0.
Please use `as_data_pronoun()` instead
This warning is displayed once per session.
Warning: `new_overscope()` is soft-deprecated as of rlang 0.2.0.
Please use `new_data_mask()` instead
This warning is displayed once per session.
Warning: The `parent` argument of `new_data_mask()` is deprecated.
The parent of the data mask is determined from either:

  * The `env` argument of `eval_tidy()`
  * Quosure environments when applicable
This warning is displayed once per session.
Warning: `overscope_clean()` is soft-deprecated as of rlang 0.2.0.
This warning is displayed once per session.
C$name = colnames(X)
  • 查看商業類別摘要
sapply(list(X=X, B=B, C=C), dim)
          X      B   C
[1,] 188593 188593 936
[2,]    936     80   4

(4) 商業類別字雲 Category Word Cloud by Businesses

  • 以下我們使用文字雲觀察商業類別之間的相似性， 使用tSNE，將X的尺度 [188593 x 936] 縮減為 [2 x 936] …
library(RColorBrewer)
library(wordcloud)
library(Rtsne)

if(!LOAD) {
  t0 = Sys.time()
  set.seed(123)
  tsneCat = Rtsne(as.matrix(t(X)), check_duplicates=F, theta=0.0, max_iter=3000)
  Sys.time() - t0   # 3.857 mins
  }
  • 在縮減尺度之中做階層式集群分析，並將同一群的商業類別以相同顏色做表示
Y = tsneCat$Y           # tSNE coordinates
d = dist(Y)             # distance matrix
hc = hclust(d)          # hi-clustering
K = 80                  # number of clusters 
C$group = g = cutree(hc,K)        # cut into K clusters
table(g) %>% as.vector %>% sort   # sizes of clusters
 [1]  2  2  4  4  4  4  4  4  5  6  6  6  7  7  7  7  7  7  8  8  8  8  8  8  8  8  8  8  9  9  9  9
[33]  9 10 10 10 10 10 10 10 11 11 11 11 11 11 11 12 12 13 13 13 14 14 14 15 15 15 15 15 15 16 16 16
[65] 16 17 17 17 17 18 18 18 20 20 21 23 23 24 28 31
  • 調整商業類別當中評論數量C$n_rev的範圍
sz = 0.7 + sqrt(C$n_rev)/500
range(sz)  
[1] 0.7211 4.5235
  • 繪製文字雲
png("fig/category.png", width=3200, height=1800)
textplot(Y[,1], Y[,2], C$name, font=2, 
         col = randomcoloR::distinctColorPalette(K)[g],
         cex = sz ) # size by no. reviews
dev.off()
png 
  2

將字雲畫在category.png裡面：

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

我們分別使用了尺度縮減和集群分析來做以上的字雲，其中 …
  ■ 尺度縮減的
    ● 原始尺度有多少個？
    ● 縮減之後剩下多少尺度？
    ● 原始尺度是什麼？換句話說，我們是根據甚麼來做尺度縮減？

  ■ 我們是根據什麼做的集群分析？
    ● 是原始尺度、還是縮減之後的尺度？
    ● 用原始和縮減尺度、會有什麼差別？


(5) 評論話題字雲 Theme Word Cloud

使用字雲觀察評論話題(Theme)之間的相似性

5.1 Average Sentiment & Empath scores per business

接下來考慮評論的話題，我們已經預先使用Stanford的Empath Text Classifier，依其預設的194種內容(Class)， 對這5,996,996篇評論分別做過評分，文集之中的每一篇評論都有194個內容評分，放在 data/Tmx.rdata 裡面；由於資料太大，我們先依商店和評論人分別對話題全做過平均。

load("data/Tmx.rdata")
  • 平均情緒評分
if(LOAD) { load("data/BCscores.rdata") } else 
{
  t0 = Sys.time()
  cat_senti = apply(X, 2, function(v) colMeans( rev[rev$bid %in% B$bid[v], 8:17] ) ) %>% t
  Sys.time() - t0
  biz_senti = aggregate(. ~ bid, data=rev[, c(2,8:17)], mean)
  Sys.time() - t0
  cat_theme = apply(X, 2, function(v) colMeans( Tmx[rev$bid %in% B$bid[v],] ) ) %>% t 
  Sys.time() - t0
  biz_theme = aggregate(as.data.frame.matrix(Tmx), list(bid = rev$bid), mean)
  Sys.time() - t0
  gc()
  save(biz_senti, biz_theme, cat_senti, cat_theme, file="data/BCscores", compress=T)
}

biz_senti - 每一個商業類別之平均情緒分數

dim(biz_senti) # biz_senti維度
[1] 188593     11

biz_theme - 每一個商業類別之平均話題主題權重

dim(biz_theme) # biz_theme維度
[1] 188593    195
5.2 話題的討論強度 - Summerize Empath Scores across Biz
themes = data.frame(name=colnames(Tmx), weight=colSums(Tmx), stringsAsFactors=F)
  • 前20個話題主題的討論強度
par(mar=c(8,4,4,2), cex=0.7)
colSums(Tmx)[1:20] %>% barplot(main="Sums of Empath Scores, Top20 Themes", las=2)

par(mar=c(2,4,4,2), cex=0.7)
themes$weight %>% barplot(main="Sums of Empath Scores")

5.4 話題文字雲 Theme Word Cloud
if(!LOAD) {
  t0 = Sys.time()
  set.seed(123)
  tsneTheme = biz_theme[,-1] %>% scale %>% as.matrix %>% t %>% 
    Rtsne(check_duplicates=F, theta=0.0, max_iter=3000)
  Sys.time() - t0  # 29.21 secs
}
Y = tsneTheme$Y         # tSNE coordinates
d = dist(Y)             # distance matrix
hc = hclust(d)          # hi-clustering
K = 40                  # number of clusters 
themes$group = g = cutree(hc,K)  # clustering for color
table(g) %>% as.vector %>% sort  # size of clusters
 [1] 2 2 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 7 7 7 7 7 8
sz = sqrt(themes$weight)/100 + 1.5
range(sz)   
[1] 1.591 5.122
  • 繪製話題文字雲
png("fig/theme.png", width=3200, height=1800)
textplot(Y[,1], Y[,2], themes$name, font=2, 
         col = randomcoloR::distinctColorPalette(K)[g],    # color by group    
         cex = sz )                                        # size by total weight
dev.off()
png 
  2

(6) 話題、類別的對應關係

6.1 TC - Theme-Category Matrix
  • 首先將評論話題與商業類別之間的關係整理成矩陣
library(d3heatmap)
#TC = apply(X, 2, function(i) 100*colMeans(E[i > 0,]) )
#dim(TC) 
#sapply(list(TC, colSums(TC), rowSums(TC)), range)
  • 使用熱圖表現出評論話題與商業類別之間的關係
library(d3heatmap)
cat_theme[1:100, 4:80] %>% t %>% d3heatmap(colors = cm.colors(13)[3:13]) 
# rev(brewer.pal(11,"Spectral"))
# 這邊取其中350個商業類別對194個話題主題來生成熱圖
# 可以用以觀察某商業類別之中，話題主題集中在哪邊
cat_theme[1:350, 1:194] %>% scale %>% d3heatmap(
  show_grid=F, 
  xaxis_font_size = "0px", xaxis_height = 10,
  yaxis_font_size = "0px", yaxis_width = 10,
  colors = brewer.pal(9,"Greens")
  )