この記事はDATUM STUDIO Advent Calendar 2016の22日目の記事である。

{LDAvis}パッケージ

LDAの結果をいい感じに可視化してくれるパッケージ。
CRANのマニュアルは以下。
https://cran.r-project.org/web/packages/LDAvis/index.html

前提条件

今回のサンプルでは日本語の形態素解析のためにRMeCabを用いています。
インストールは下記参照。
http://rmecab.jp/wiki/index.php?RMeCab

動作確認

Mac環境での動作を確認しています。
```

sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)
```

解析に用いた文章

例として、3つのニュース記事のタイトルと本文を処理しています。
```
$ cat text_1.txt
博多駅前再び道路７センチほど沈み込む通行止めに
２６日未明、福岡市のＪＲ博多駅前の大規模に道路が陥没した現場付近で、再び道路が最大で深さ７センチほど沈んでいるのが見つかり、警察は周辺の交通を規制して、詳しい状態などを調べています。
・・・

$ cat text_2.txt
陥没めど立たぬ休業補償　飲食店など博多駅前事業者　福岡市に問い合わせ５０件

福岡市のＪＲ博多駅前の道路陥没事故で休業を余儀なくされた飲食店などの事業者から損失補償に関する問い合わせが市に相次ぎ、１１日までに５０件を超えた。
・・・

$ cat text_3.txt
陥没周辺の建物は？　福岡市「倒壊の恐れなし」　専門家「地震の揺れに注意を」

陥没事故の周辺では、建物の倒壊など二次被害も懸念された。福岡市は８日、陥没箇所を中心に東西約４００メートル、南北約１５０メートルにある４２棟で応急危険度判定を実施。「地上部や地中のくいに傾きが見られず、現状で倒壊の恐れがある建物はない」とする調査結果を公表した。
・・・
```

実行例

以下に実行例を示します。

日本語テキストの形態素解析

library(RMeCab)

## 3つの文章からそれぞれ名詞、形容詞のみ残す
## 文章その1
tmp_doc.1 <- NULL
tmp_doc.1 <- RMeCabText("text_1.txt")
doc.1 <- NULL
for (i in 1:length(tmp_doc.1)) {
    if (tmp_doc.1[[i]][2] %in% c("名詞", "形容詞")) {
        doc.1 <- c(doc.1, paste(tmp_doc.1[[i]][1], sep = "", collapse = " "))
    }
}

## 文章その2
tmp_doc.2 <- RMeCabText("text_2.txt")
doc.2 <- NULL
for (i in 1:length(tmp_doc.2)) {
    if (tmp_doc.2[[i]][2] %in% c("名詞", "形容詞")) {
        doc.2 <- c(doc.2, paste(tmp_doc.2[[i]][1], sep = "", collapse = " "))
    }
}

## 文章その3
tmp_doc.3 <- RMeCabText("text_3.txt")
doc.3 <- NULL
for (i in 1:length(tmp_doc.3)) {
    if (tmp_doc.3[[i]][2] %in% c("名詞", "形容詞")) {
        doc.3 <- c(doc.3, paste(tmp_doc.3[[i]][1], sep = "", collapse = " "))
    }
}

#各文章をリスト化
doc.list <- NULL
doc.list <- list(doc.1, doc.2, doc.3)
names(doc.list) <- c("doc1", "doc2", "doc3")

ldaの実行

library(lda)
# ターム行列作成
term.table <- table(unlist(doc.list))
term.table <- sort(term.table, decreasing = TRUE)

# 単語一覧
vocab <- NULL
vocab <- names(term.table)

# 単語のインデックス化とそれぞれの出現回数をカウント
get.terms <- function(x) {
  index <- match(x, vocab)
  index <- index[!is.na(index)]
  rbind(as.integer(index - 1), as.integer(rep(1, length(index))))
}
documents <- NULL
documents <- lapply(doc.list, get.terms)

## パラメータ設定
D <- length(documents)  # ドキュメント数
W <- length(vocab)  # 単語の種類
doc.length <- sapply(documents, function(x) sum(x[2, ]))  # それぞれの文章の単語数
N <- sum(doc.length)  # 総単語数
term.frequency <- as.integer(term.table)  # 頻度行列


# 学習用パラメータ
K <- 10 ## トピック数
G <- 200 ## 反復回数
alpha <- 0.02 
eta   <- 0.02

## ldaの実行
fit <- lda.collapsed.gibbs.sampler(documents = documents, 
                                   K = K,
                                   vocab = vocab, 
                                   num.iterations = 25,  # 繰り返し数
                                   alpha = alpha, # ディリクレ過程のハイパーパラメータα
                                      eta = eta, # ディリクレ過程のハイパーパラメータη
                                      compute.log.likelihood=TRUE)

# 各トピックのトップ3ワード
top.words <- top.topic.words(fit$topics, 3, by.score = TRUE)
print(top.words)

LDAvisで可視化(本題)

library(LDAvis)

# パラメータ
display_word_num = 20 #表示する単語数
theta      <- t(apply(fit$document_sums + alpha, 2, function(x) x/sum(x)))
phi        <- t(apply(t(fit$topics) + eta, 2, function(x) x/sum(x)))
doc.length <- sapply(documents, function(x) length(x)) 

# 
json <- createJSON(phi            = phi, 
                   theta          = theta, 
                   doc.length     = doc.length, 
                   vocab          = vocab, 
                   term.frequency = term.frequency,
                   R = display_word_num
                   )

## 可視化。古いvisフォルダがある場合は削除して実行
serVis(json, out.dir = 'vis', open.browser = T )

出力結果

インタラクティブにUIの動作が確認可能。
道路、原因、保証などいくつかのトピックが確認できる。

マウスオーバーでそれぞれのトピックの中身を確認できる。

なおserVisで指定したout.dirにLDAvisの結果が一式出力されます。
Firefoxなどのブラウザでindex.htmlを開けばRのない環境でも結果の確認が可能になります。

最後に

来月から福岡で働くことになるのだが、飲みすぎて穴に落ちないことをただ祈るばかりである。

enjoy!