【Rで自然言語処理】トピックモデルの階層構造をどうにか可視化したい。
もう7月ですね。
統計検定、終わっちゃいましたね。
結果云々を一切書いてないってことはお察しということで。
※というか、6月末だと思い込んでいて、余裕ぶっこいていたら申し込みすらすっかり忘れていた、というオチです。
はい、というわけで、今回はトピックモデルのお話です。
ふと思ったらトピックモデルを本来の用途で使った記事を書いたことがなかったなぁ、と。
wanko-sato.hatenablog.com
こんなのは書いたけれども、扱ったのは数字の集合だったわけで、さすがにこれで「トピックモデルやりました」とは言いづらいかなぁ、と。
とはいえ、トピックモデルの詳しい話はしません。理論的な背景やら実際の計算方法やらは自分なんかが説明するよりも書籍を買って読んでいただく方がずっと詳しいです。
- 作者: 岩田具治
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (2件) を見る
トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)
- 作者: 佐藤一誠,奥村学
- 出版社/メーカー: コロナ社
- 発売日: 2015/03/13
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (5件) を見る
今回、この記事を書こうと思ったモチベーションは、
- LDAはなかなか性能が良いし、なかなかきれいに分類してくれてうれしい
- でも、トピックモデルの解釈は意外に面倒
- HDP-LDAでトピック数の自動決定は可能だけれど、それもハイパーパラメータに依存するし、どうにもいまいちな感じがしている
- できれば、大分類＞中分類＞小分類くらいな感じの階層構造で分けたい
- それを可視化できれば、もっと結果の解釈が楽になるのでは？
といったところです。ここのところいろいろ考えているのが分析結果の可視化でして。ぱっと見て概要がすぐわかり、ドリルダウンも簡単、細かい情報は元のデータのリストを参照してね、というふうにしたいのです。分析結果は分析した人がわかればいいのではない、いろんな人がいろんな角度から考察できるようにアシストできる形で示すべきだ、と考えるからです。
能書きはここまでにして、さっそく中身に入りましょう。
- 今回やること
- PubMedの検索結果をダウンロードする
- MEDLINEのデータを読み込む
- テキスト情報を加工し、LDAでモデリングする
- LDAの結果を可視化する（１）
- LDAの結果を可視化する（２）
- javascriptで可視化する
- 今後
今回やること
今回は次のようなことをやります。
PubMedの検索結果をダウンロードする
PubMedはバイオメディカル系の文献データベースです。論文のタイトル、掲載ジャーナル、出版年月日、著者、要旨などのデータをMEDLINEと呼ばれる構造化されたテキストデータとしてダウンロードできます。今回はそのMEDLINEのデータから論文の要旨等を利用して論文の分類をトピックモデルで行おう、というのが基本的な方針です。
ちなみに、XML形式でも出力できるのですが、階層構造が複雑で扱いが面倒なので、個人的にはMEDLINE形式がおすすめです。
今回は以下のキーワードで検索した論文のうち、出版日が2017年6月1日以降のものを使用しました。
- allergy
- breast cancer
- colorectal cancer
- gastric cancer
- psychiatry
アレルギー、がん、精神疾患と大きく分け、がんを乳がん、大腸がん、胃がんの3種類に分類しました。この階層構造がうまく表現できるか、というのが今回の主眼になります。
ちなみに検索結果をダウンロードするには、検索結果ページの右上にある「send to」をクリック、「File」を選択し、Formatを「MEDLINE」にして「Create File」をクリックすればOKです。
ダウンロードしたファイルはテキストファイルです。中身を見ていただくと、タグとそのタグの詳細内容、という単純な形式になっています。タグの細かい説明はWikiやPubMedのヘルプをご覧ください。
では、次に行きましょう。
MEDLINEのデータを読み込む
Text Mining PubMed: Parsing Medline Files in Rmlbernauer.wordpress.com
データ読み込みの部分はほぼ丸パクリさせていただきました。
MELINE読み込み用の関数を定義する
## define function
medline <- function(file_name){
lines <- readLines(file_name)
medline_records <- list()
key <- 0
record <- 0
for(line in lines){
header <- sub(" {1,20}", "", substring(line, 1, 4))
value <- sub("^.{6}", "", line)
if(header == "" & value == ""){
next
}
else if(header == "PMID"){
record = record + 1
medline_records[[record]] <- list()
medline_records[[record]][header] <- value
}
else if(header == "" & value != ""){
medline_records[[record]][key] <- paste(medline_records[[record]][key], value)
}
else{
key <- header
if(is.null(medline_records[[record]][key][[1]])){
medline_records[[record]][key] <- value
}
else {
medline_records[[record]][key] <- paste(medline_records[[record]][key], value, sep=";")
}
}
}
return(medline_records)
}この関数を定義し、実行すればOK。
複数のファイルを読み込む
今回は5つのファイルをいっぺんに読み込むため、次のようにしました。
## read Medline
fileVec <- list.files(file.path(getwd(),"R","topicmodel"),
pattern=".txt",
full.names = T)
categoryVec <- list.files(file.path(getwd(),"R","topicmodel"),
pattern=".txt")
dataMed <- lapply(fileVec,medline)まずファイルのフルパスを取得してベクトルに格納します。
次のcategoryVecはどの検索ワードで引っ張ってきた論文かを確認するための識別用の情報で、後で使います。
最後にlapplyでmedline関数を呼び出し、Listにします。
なお、medline関数は読み込んだファイルをList形式で格納するので、ListのListになっています。そのままでも良いのですが、ちょっと扱いづらいので使うデータだけをデータフレームの形に直します。
dataMedList <- lapply(seq(1,length(dataMed)),function(x){
res <- lapply(seq(1,length(dataMed[[x]])),function(y){
out <- data.frame(title=dataMed[[x]][[y]]$TI,
abst=dataMed[[x]][[y]]$AB,
category=categoryVec[x])
return(out)
})
return(do.call(rbind,res))
})
dataMedDF <- do.call(rbind,dataMedList)これで、論文のタイトルと要旨、ならびにカテゴリを格納したデータフレームができました。
※このデータフレームの出番は後半です。LDAのモデリングまではほとんど使いません。
付加情報をもってくる
PubMedに掲載されている論文には付加情報としてMeshというキーワードのようなものが割り当てられています。論文につけるキーワードは著者が任意で決められるため、同じ意味でも異なる表現になることが多々あります。それでは分類に不向きなので、PubMedが設定した統制用語を割り当てています。その統制用語がMeshです。
Meshにも階層構造があって、疾患領域であったり薬剤の分類であったり、使おうと思えばかなり有用な情報にできるのですが、今回はそこまではやらず、あくまで割り当てられたMeshを要旨にくっつけるだけにとどめます。
dataMH <- sapply(dataMed,function(i){
sapply(i,function(x){
if(all(!(names(x) %in% "MH"))){
strMH <- c()
}else{
strMH <- strsplit(x$MH,";")
}
strMHgsub <- gsub(" ","_",strMH[[1]])
strMHpaste <- paste(strMHgsub,collapse = " ")
return(strMHpaste)
})
})ここで注意！MEDLINEのタグには"MH"と"MHDA"という二つのタグがあります。MHはMeshそのものなのですが、MHDAはMeshが付与された日付を格納しています。で、Rの仕様なのか、リスト要素をnameで指定する場合にMHを使うとMHDAまで引っ張ってきてしまうのです。そういうわけで、上のコードのようにMHに完全一致するnameをもっている場合にのみ、MHを取得する、という処理をしているのです。ちょっと面倒ですが、ノイズを除くために必要な処理です。
テキスト情報を加工し、LDAでモデリングする
モデリングまでのコードを一気に載せます。
tokenize、vocablizeあたりは"text2vec"定型処理なので、そのままコピペしてもらって大丈夫です。
library(textmineR)
library(text2vec)
library(tm)
library(topicmodels)
# create vector of abstracts
allTexts <- sapply(dataMed,function(i){
sapply(i,function(x){
ti <- gsub("\\[|\\]","",x$TI)
if(all(!(names(x) %in% "MH"))){
strMH <- c()
}else{
strMH <- strsplit(x$MH,";")
}
strMHgsub <- gsub(" ","_",strMH[[1]])
strMHpaste <- paste(strMHgsub,collapse = " ")
paste(ti,x$AB,strMHpaste,collapse=" ")
})
})
allTexts <- unlist(allTexts)
# preprocess
sw <- c("i", "me", "my", "myself", "we", "our", "ours",
"ourselves", "you", "your", "yours", tm:::stopwords("English"))
preText <- tolower(allTexts)
preText <- tm::removePunctuation(preText)
preText <- tm::removeWords(preText,sw)
preText <- tm::removeNumbers(preText)
preText <- tm::stemDocument(preText, language = "english")
# tokenize(split into single words)
it <- itoken(preText,
preprocess_function = tolower,
tokenizer = word_tokenizer,
ids = abstID)
# delete stopwords and build vocablary
vocab <- create_vocabulary(it,
stopwords = sw)
# word vectorize
vectorizer <- vocab_vectorizer(vocab)
# create DTM
dtm <- create_dtm(it, vectorizer)
# Term frecency
TDF <- TermDocFreq(dtm)
# modeling by package "topicmodels"
model <- LDA(dtm, control=list(seed=37464847),k = 20, method = "Gibbs")
termsDF <- get_terms(model,100)
topicProbability <- data.frame(model@gamma,dataMedDF$category)今回はパッケージ"topicmodels"のLDAを使っています。text2vecにもLDAは入っているのですが、各文書のトピックへの割り当て確率thetaが直接結果に含まれておらず、transformしないと出せないため、ちょっと面倒でやめにしました。そこで、どうせなら、とLDA以外の種々のモデリング手法もカバーしているtopicmodelsパッケージを使うことにしたのです。
今回はトピック数kを20で固定しました。一般にトピック数が少ないと分類が粗くなり、トピック数が多いと分類が細かくなります。といって、だんだんトピック数を多くしていけば階層構造になるか、というとそんなこともありません。あくまで所与のトピック数にあった確率分布を構築するためです。
コードの最後の方でtopicProbabilityを出力しています。これは各文書が所与のトピック20個のうち、どのトピックにどれくらいの確率で割り当てられるか、を示したデータです。
> head(topicProbability)
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
1 0.020454545 0.02045455 0.156818182 0.011363636 0.01136364 0.01136364 0.034090909 0.01590909 0.03409091 0.15681818
2 0.017156863 0.01715686 0.012254902 0.026960784 0.04656863 0.01715686 0.017156863 0.02205882 0.03676471 0.02696078
3 0.029296875 0.02539062 0.033203125 0.056640625 0.04101562 0.01757812 0.009765625 0.01367188 0.08007812 0.01367188
4 0.008650519 0.01211073 0.008650519 0.008650519 0.01211073 0.01903114 0.008650519 0.01557093 0.01211073 0.01903114
5 0.014204545 0.02556818 0.019886364 0.019886364 0.02556818 0.03693182 0.019886364 0.01420455 0.01988636 0.05965909
6 0.025193798 0.06395349 0.017441860 0.021317829 0.06395349 0.01356589 0.029069767 0.03682171 0.05232558 0.02906977
X11 X12 X13 X14 X15 X16 X17 X18 X19 X20
1 0.01136364 0.02954545 0.020454545 0.29772727 0.05681818 0.015909091 0.02045455 0.03409091 0.01136364 0.02954545
2 0.05637255 0.01225490 0.031862745 0.50735294 0.05637255 0.012254902 0.04166667 0.01225490 0.01715686 0.01225490
3 0.32617188 0.09570312 0.033203125 0.08398438 0.02148438 0.037109375 0.01757812 0.02929688 0.02148438 0.01367188
4 0.61764706 0.02249135 0.008650519 0.01557093 0.12283737 0.008650519 0.01903114 0.02941176 0.01903114 0.01211073
5 0.01420455 0.02556818 0.014204545 0.07670455 0.12215909 0.014204545 0.01420455 0.03693182 0.03693182 0.38920455
6 0.01744186 0.02131783 0.040697674 0.02519380 0.42441860 0.017441860 0.02131783 0.01356589 0.05232558 0.01356589
dataMedDF.category
1 pubmed_result_allergy.txt
2 pubmed_result_allergy.txt
3 pubmed_result_allergy.txt
4 pubmed_result_allergy.txt
5 pubmed_result_allergy.txt
6 pubmed_result_allergy.txtこの中でいうと、文書1はのトピック14への割り当てが30%程度となっており、最も高い値となっている、という具合。
もし、文書同士が似ているとするならば、この各トピックへの割り当て具合も似てくるはずです。そこで、まずはこのデータを使って類似度を出してみよう、という発想になるわけです。
LDAの結果を可視化する（１）
というわけで、ぱっと思いつくところで、topicProbabilityからユークリッド距離を求めて可視化してみましょう。
ユークリッド距離を求める
これは簡単です。
# similarity edist <- dist(model@gamma,method="euclidean") edistMat <- as.matrix(edist) point <- cmdscale(edistMat) plot(point)
2次元平面上にプロットする
色分けは次のようになっています。
- 赤 : allergy
- 緑 : breast_cancer
- 青 : colorectal_cancer
- 水 : gastric_cancer
- 紫 : psychiatry
うーん、紫（というかピンク）がごしゃっと重なっていてわかりづらい絵になってしますねぇ。。。
3次元表示にするともうちょいましでしょうか。
3次元空間上にプロットする
library(rgl) point <- cmdscale(edistMat,3) plot3d(point)
これでだいぶわかりやすくなりました。
ピンクのpsychiatryが左下の方に分離されています。左上方向に青のcolorectal cancer、右下方向に緑のbreast cancerが分離され、水色のgastric cancerは青に近い位置に配置されています。赤のallergyは角の真ん中らへんに集中している感じでしょうか。
ここから、psychiatryはほかの論文とはっきり分離されていることがわかります。また、各がんはpsychiatryと異なる分類になっており、さらに消化器系と乳がんが分けられている感じです。
こうしてみると、LDAでずいぶんときれいに文書の分類ができるものだなぁ、と感心してしまいます。
なのですが。
- じゃあ、この角の先端と中心とはどう違うのだろうか？
- psychiatryとそれ以外は大きく二つの方向に分かれているが、それはどのような分類になっているのだろうか？
といったことはここからはぱっと読み取れません。
そういう理由もあり、トピックモデルの結果を階層構造にしてはどうか、という発想に至ったわけなのです。
LDAの結果を可視化する（２）
大分類＞中分類＞小分類という階層構造をどう表現しようか？と考えたときに見つけたのがこれです。
Zoomable Circle Packing - bl.ocks.org
各文書がどのクラスターに属するか、を階層構造で表現できるもの。javascriptでの実装になりますが、せっかくだし、と思い頑張ってやってみました。
このグラフを表示するためには、json形式の階層構造のデータが必要になります。そこで、hclustで階層クラスタリングしたデータをjson形式に変換していく一連の流れをやってみました。
※細かいところまではぜんぜんできてません。
階層クラスタリング
実験もかねて、とりあえず手っ取り早いhclustをやります。
# hierlacy cluster hclustRes <- hclust(edist) plot(hclustRes)
数が多いとわけわかんないです。
これをもっとわかりやすくするのが目標です。
階層構造のデータを作る
ここからが今回の記事のポイントです。
このあたりを参考にしながらやっていきます。たぶんもっといいやり方があると思うんですが、とりあえず今回はこれで。
使うライブラリは以下の通りです。
library(RJSONIO) #RJSONIOがポイント
library(treemap)
library(d3treeR) # devtools::install_github("timelyportfolio/d3treeR")では、いきましょう。
hierStr <- cutree(hclustRes,h=c(0.75,0.5,0.25,0))
まず、デンドログラムを適当な深さで切っていきます。今回は0.75,0.5,0.25としました。きちんと分割する深さを評価してやればいいのでしょうけれど、今回は実験ということで。
ちなみに0は一番下、分割前の全要素になります。
hierStr <- as.data.frame(cbind(rep("root2",nrow(edistMat)),hierStr))
colnames(hierStr) <- c("root","level1","level2","level3","level4")頭に大本の根っこになる"root2"をつけました。
こんな感じのデータフレームになります。
> head(hierStr) root level1 level2 level3 level4 1 root2 1 1 1 1 2 root2 1 1 2 2 3 root2 2 2 3 3 4 root2 2 2 4 4 5 root2 3 3 5 5 6 root2 4 4 6 6
各要素がどの枝に属しているか、が深さ毎にわかるようなデータです。
これに付加情報をつけてざくざくやっていくわけです。
ツリーマップに変換する
最初に付加情報をくっつけておきます。
※実は今回は使えていませんが。。。
hierStr <- data.frame(hierStr,
size=c(1000),
id=hclustRes$order,
name=dataMedDF[hclustRes$order,]$title,
category=dataMedDF[hclustRes$order,]$category)文書ID、タイトル、カテゴリをくっつけています。sizeは表示するときのノードの大きさです。今回はすべて同じ値にしています。
※たとえばこれを論文掲載誌のインパクトファクターにする、とかいろいろ工夫できるはずです。
次に、このデータフレームをツリーマップに変換します。
なんでいちいちツリーマップに変換するのか？というのも、このやり方が一番手っ取り早かったからです。もとのデータフレームをあれこれいじくってjsonにしようと思ったのですが、意外に手数がかかるため、だったらなるべく少ない手数でできる方法でいこう、という方針で至った結論です。
indexList <- c("root","level1","level2","level3","level4")
treedat <- treemap(hierStr, index=indexList, vSize='size')
treedat <- treedat$tm
はい、こんな感じのツリーマップができました。
これ見ていただくと、たしかにやろうと思っている階層構造ができているわけで。
> head(treedat)
$tm
root level1 level2 level3 level4 vSize vColor stdErr vColorValue level x0 y0 w
1 root2 1 1 1 1 1000 1 1000 NA 5 0.85561464 0.06734416 0.03919965
2 root2 1 1 1 <NA> 1000 1 1000 NA 4 0.85561464 0.06734416 0.03919965
3 root2 1 1 2 18 1000 1 1000 NA 5 0.78505527 0.16991928 0.02351979
4 root2 1 1 2 19 1000 1 1000 NA 5 0.78505527 0.11863172 0.02351979
5 root2 1 1 2 2 1000 1 1000 NA 5 0.80857506 0.16991928 0.02351979
6 root2 1 1 2 36 1000 1 1000 NA 5 0.83209485 0.16991928 0.02351979
7 root2 1 1 2 39 1000 1 1000 NA 5 0.80857506 0.11863172 0.02351979こういうデータです。"root2"をトップにして各ノードがどの階層のどのクラスターに属しているか、がデータになっています。で、このtreemapのデータを一発でjsonに変換する関数があるわけです。
json形式に変換する
res <- d3treeR:::convert_treemap(treedat, rootname="flare") json <- RJSONIO::toJSON(res, auto_unbox = TRUE) writeLines(json,file.path(getwd(),"R","topicmodel","treeJson.json"))
d3treeRパッケージにconvert_treemapという関数がありまして、ツリーマップのデータを階層構造のList形式に変換してくれるのです。で、RJSONIOでjson形式に変換、最後にファイル出力して完成、です。
ちなみに、出力されたjsonのファイルは
{
"name": "flare",
"id": 1,
"size": null,
"children": [
{
"name": "root2",
"color": "#00BCD2",
"h": 1,
"id": 2,
"size": 8.29e+05,
"stdErr": 8.29e+05,
"vColor": 829,
"vSize": 8.29e+05,
"w": 1,
"x0": 0,
"y0": 0,
"children": [
{
"name": "1",
"color": "#54A249",
"h": 0.23283,
"id": 3,
"size": 24000,
"stdErr": 24000,
"vColor": 24,
"vSize": 24000,
"w": 0.12434,
"x0": 0.69432,
"y0": 0,
"children": [
{
・・・こんな感じです。インデントもちゃんとしてないし、余計な要素もいっぱいくっついていますが、とりあえず表示できればよし！の精神で次に進みたいと思います。
javascriptで可視化する
javascriptのコードはほぼ丸パクリです。
<!DOCTYPE html> <meta charset="utf-8"> <style> .node { cursor: pointer; } .node:hover { stroke: #000; stroke-width: 1.5px; } .node--leaf { fill: white; } .label { font: 11px "Helvetica Neue", Helvetica, Arial, sans-serif; text-anchor: middle; text-shadow: 0 1px 0 #fff, 1px 0 0 #fff, -1px 0 0 #fff, 0 -1px 0 #fff; } .label, .node--root, .node--leaf { pointer-events: none; } </style> <svg width="750" height="750"></svg> <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.9.1/d3.min.js"></script> <script> var svg = d3.select("svg"), margin = 20, diameter = +svg.attr("width"), g = svg.append("g").attr("transform", "translate(" + diameter / 2 + "," + diameter / 2 + ")"); var color = d3.scaleLinear() .domain([-1, 5]) .range(["hsl(152,80%,80%)", "hsl(228,30%,40%)"]) .interpolate(d3.interpolateHcl); var pack = d3.pack() .size([diameter - margin, diameter - margin]) .padding(2); var select = document.querySelector('#selectFileInputFile'); d3.json("treeJson.json", function(error, root) { if (error) throw error; root = d3.hierarchy(root) .sum(function(d) { return d.size; }) .sort(function(a, b) { return b.value - a.value; }); var focus = root, nodes = pack(root).descendants(), view; var circle = g.selectAll("circle") .data(nodes) .enter().append("circle") .attr("class", function(d) { return d.parent ? d.children ? "node" : "node node--leaf" : "node node--root"; }) .style("fill", function(d) { return d.children ? color(d.depth) : null; }) .on("click", function(d) { if (focus !== d) zoom(d), d3.event.stopPropagation(); }); var text = g.selectAll("text") .data(nodes) .enter().append("text") .attr("class", "label") .style("fill-opacity", function(d) { return d.parent === root ? 1 : 0; }) .style("display", function(d) { return d.parent === root ? "inline" : "none"; }) .text(function(d) { return d.data.name; }); var node = g.selectAll("circle,text"); svg .style("background", color(-1)) .on("click", function() { zoom(root); }); zoomTo([root.x, root.y, root.r * 2 + margin]); function zoom(d) { var focus0 = focus; focus = d; var transition = d3.transition() .duration(d3.event.altKey ? 7500 : 750) .tween("zoom", function(d) { var i = d3.interpolateZoom(view, [focus.x, focus.y, focus.r * 2 + margin]); return function(t) { zoomTo(i(t)); }; }); transition.selectAll("text") .filter(function(d) { return d.parent === focus || this.style.display === "inline"; }) .style("fill-opacity", function(d) { return d.parent === focus ? 1 : 0; }) .on("start", function(d) { if (d.parent === focus) this.style.display = "inline"; }) .on("end", function(d) { if (d.parent !== focus) this.style.display = "none"; }); } function zoomTo(v) { var k = diameter / v[2]; view = v; node.attr("transform", function(d) { return "translate(" + (d.x - v[0]) * k + "," + (d.y - v[1]) * k + ")"; }); circle.attr("r", function(d) { return d.r * k; }); } }); </script>
これはこのままでは表示できません。というのも、javascriptでローカルの外部ファイルを読み込む場合、ただindex.htmlをブラウザに表示させるだけではだめだからです。なので、今回はpythonから簡易サーバを立てる方法をとりますた。参考は以下のサイトです。
※pythonがインストールされていないとできません。入っていない場合はインストールしてください。
python -m http.server 8888
で簡易サーバを立ち上げ、ブラウザのアドレスバーに"localhost:8888"と入力、目的のindex.htmlが保存されているディレクトリまで移動していくと・・・
できたーーー！
ちゃんとズームもできます！！！
いやー、javascriptってすごいですね！
※コードは丸パクリだけどな。
※はてなブログにjavascriptが埋め込める、という噂を聞いたことがあるのですが、今回はそこまでやりません。
今後
とりあえず当初の目的は達成できたわけなんですが、このまんまだとやっぱりよくわからないままなので、付加情報もちゃんと表示できるようにしたい。
また、階層クラスタリングが本当に妥当なのか、分割の仕方が本当に妥当なのか、をちゃんと検証する必要がある。
などなど、いろいろ課題もあるわけです。
が、だいぶ長くなってしまったのでひとまずこれでお開き。
