【Rで自然言語処理】トピックモデルによる文書分類結果の妥当性を検証する。

前回、トピックモデルのLDAを用いて文書を分類し、その結果をjavascriptで可視化する、ということをやりました。
wanko-sato.hatenablog.com
ただ、その結果が本当に妥当なのか、正直自信がありませんでした。一応は3Dプロットで、検索ワードと文書分類がほどほどかみ合っていることは確認しましたが、あまりにざっくりしすぎてきちんとした検証を行っていませんでした。
というわけで、前回使ったデータと結果をもとに、得られた結果が実際のところどんなものなのか、じっくり見ていきたいと思います。

前回のおさらい

前回やったことは大体こんなところです。

1. PubMedから異なる検索ワード5種類で引っ張ってきた論文のabstractをLDAにかける
2. LDAのトピック数を20に固定したときに得られた各トピックへの割り当て確率をもとに文書間のユークリッド距離を算出する
3. 距離行列を3次元座標に落とし込み、可視化する
4. 距離行列を使って階層的クラスタリングを行う
5. デンドログラムを任意の数か所で切断し、クラスターの階層構造をツリーマップに落とし込む
6. ツリーマップをjson形式のデータに吐き出し、javascriptで可視化する

なぜこんなことをしようと思ったのか、そもそものモチベーションは「トピックモデルのLDAは非常に良い手法だと思うのだけれど、結果の解釈がいかんせん煩雑なため、ビジュアル的かつインタラクティブに結果を眺められないだろうか？」というところでした。

3番目、距離行列を3次元座標に落とし込み（Rのcmdscale関数を使用）プロットしたものが上の図です。色の分け方は

• 赤 : allergy
• 緑 : breast_cancer
• 青 : colorectal_cancer
• 水 : gastric_cancer
• 紫 : psychiatry

となっております。大まかにみて、左下に伸びている紫のpsychiatryと右方向に伸びている3種のcancer、あちこちに混ざっているallergy、という構図が見て取れます。こいつを階層的クラスタリングにかけたとき、はたしてうまいこと分類されているのだろうか？というのが、今回検証したいことの一番目です。

以下の記事を読み進めていくにあたり、大前提として前回のRコードを一通り実行し、本記事中にいきなり登場するデータフレームやリストがすでに生成されているものとします。

階層的クラスタリング

コードの修正をいくつか

まず、前回のコードで

# hierlacy cluster
hclustRes <- hclust(edist)
plot(hclustRes)

としておりましたが、本来は

# hierlacy cluster
hclustRes <- hclust(edist,method="ward.D2")
plot(hclustRes)

とすべきでした。hclustはデフォルトでmethod="complete"となっていますが、できればウォード法を使いたいところです。ウォード法でhclustした場合、デンドログラムの高さが必ずしも1にはならないため、デンドログラムを任意の箇所で切断する場合、

cutH <- c((max(hclustRes$height)/4)*3,
          (max(hclustRes$height)/4)*2,
          (max(hclustRes$height)/4)*1,
          0)
hierStr <- cutree(hclustRes,h=cutH)
hierStr <- as.data.frame(cbind(rep("root2",nrow(edist)),hierStr))
colnames(hierStr) <- c("root","level1","level2","level3","level4")

とする必要があります。上記のコードではとりあえず全体を4分割しています。もっと良い切断の仕方があるとは思いますが、ひとまずはこれで。

hclustの出力結果に関する注意点

あれこれ実験していて「？」と思うことが多数あり、よくよく検討してみたらわかったこと。
上のコードでhierStrというデータフレームにhclustの結果のうち、各切断点でクラスタリングした際のクラスターの割り当てを格納しています。

> head(hierStr)
   root level1 level2 level3 level4 size id
1 root2      1      1      1      1 1000  1
2 root2      1      1      1      2 1000  2
3 root2      2      2      2      3 1000  3
4 root2      2      2      2      4 1000  4
5 root2      1      3      3      5 1000  5
6 root2      2      4      4      6 1000  6

中身を見るとなんの変哲もないデータフレームだと、すっかり思い込んでいました。

> is.numeric(hierStr$level1)
[1] FALSE

！！！
格納されているのが数値かと思ったら、実は違っていて、

> is.factor(hierStr$level1)
[1] TRUE

こいつらfactorなんですよ。
まぁ、level1～level3まではfactorで表示されている値と実際の中身が合致しているようなので影響ないのですが、最大の問題はlevel4でして。
この場合のlevel4はデンドログラムの一番下、つまり各文書の番号に相当します。最終的なjsonへの出力でいえば各ノードに相当します。で、何が問題かというと、

> head(hierStr$level4)
[1] 1 2 3 4 5 6
829 Levels: 1 10 100 101 102 103 104 105 106 107 108 109 11 110 111 112 113 114 115 116 117 118 119 12 120 ... 99
> head(as.numeric(hierStr$level4))
[1]   1 112 223 334 445 556

となっており、factorの表示上の値とそこに実際に格納されている値が異なっているのです。このままだと、hierStr$level4をそのまま使って文書番号と紐づけしようと思うとおかしなことになってしまうのです。実際には格納されている値ではなく、factorの表示上の値が文書番号と一致しているので、

as.numeric(as.character(hierStr$level4))

という処理が必要になってくるのです。
最初これになかなか気づかず、おかしな結果になって混乱していたのでした。

javascriptでの可視化をカラフルに

前回示した結果は

これですが、すべてのノードが白になっており、なんのこっちゃわからない状態でした。
そこで、各ノードに3Dプロットでやったような色付けをやりたいと思ったわけです。

stackoverflow.com

この質問に対する回答がそのまんまで、javascriptの方をちょろっといじってあげれば良いわけです。
加えて、読み込ませるjsonの方にも色情報を載せてあげれば良いわけで、

# set category color
colSet <- c("#FF0000","#00FF00","#0000FF","#33FFFF","#FF33FF","#FFFF33","FFCCCC")
catNum <-as.numeric(dataMedDF$category)
catCol <- colSet[catNum]

hierStr <- data.frame(hierStr,
                      size=c(1000),
                      id=as.numeric(as.character(hierStr$level4)))

# write json
library(treemap)
indexList <- c("root","level1","level2","level3","level4")
treedat <- treemap(hierStr, index=indexList, vSize='size')
treedat <- treedat$tm
treedat$color <- catCol[as.numeric(as.character(treedat$level4))]

res <- d3treeR:::convert_treemap(treedat, rootname="flare")
json <- RJSONIO::toJSON(res, auto_unbox = TRUE)
json <- gsub("color","fill",json)
writeLines(json,file.path(getwd(),"R","topicmodel","treeJson.json"))

とします。ツリーマップのデータのcolor要素に各カテゴリの色を設定し、json上でcolorとなっているところをfillに置換するだけです。これで、各ノードに各カテゴリに対応する色を割り当てることができました。

では、見てみましょう。

おお～～。
うまくいきました。
色分けは

• 赤 : allergy
• 緑 : breast_cancer
• 青 : colorectal_cancer
• 水 : gastric_cancer
• 紫 : psychiatry

となっています。psychiatryがしっかり分離されていること、breast_cancer単独のクラスタが真ん中上にできていること、3種のcancerが混ざったクラスタがいくつかあること、などなど、あれやこれや混ざってはいるものの、それなりに各カテゴリにうまいこと分離できていることがわかります。自分でやってて軽く感動しました。

検証

さて、ここまででひとまず、LDAによる各トピックへの割り当て確率を使った階層的クラスタリングはなかなか良いクラスタリングをしてくれることがわかりました。じゃあ、この分類された各クラスターに含まれている文書集合は一体何ものなのか？という当然の疑問が出てきます。LDAはどうも、そのあたりの解釈が面倒なのです。

各クラスターに割り当てられたトピックを探す

前回も示しましたが、各トピックへの割り当ては

> 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

このような行列になっています。正直、ここからどの文書がどのトピックに割り当てられているかを探すのはかなりしんどいです。なので、ぱっと見やすいグラフにしたいのが人情です。というわけで、こんな風にしてみてはどうでしょうか。

probL1 <- as.data.frame(sapply(seq(1,length(levels(hierStr$level1))),function(x){colMeans(topicProbability[hierStr$level1==x,1:20])}))
probL2 <- as.data.frame(sapply(seq(1,length(levels(hierStr$level2))),function(x){colMeans(topicProbability[hierStr$level2==x,1:20])}))
probL3 <- as.data.frame(sapply(seq(1,length(levels(hierStr$level3))),function(x){colMeans(topicProbability[hierStr$level3==x,1:20])}))

各level毎に各トピックへの割り当ての平均を算出し、箱ひげ図を描いて各クラスタにおける外れ値にあたるトピックを抽出するのです。

# check topic
bp <- boxplot(probL1)
outTopic <- data.frame(x=bp$group,
                       y=bp$out)
whTopic <- sapply(seq(1,nrow(outTopic)),function(x){
  which(outTopic[x,2]==as.vector(probL1[,outTopic[x,1]]))
})
outTopic <- data.frame(outTopic,
                       no = whTopic,
                       txt = paste("topic",whTopic,sep=""))
boxplot(probL1)
text(outTopic$x,outTopic$y,outTopic$txt)

boxplot関数で吐き出したbpには外れ値がどれか、の情報が含まれているので、外れ値がどのトピックに当たるのかをtextで表示してやればよろしいわけです。上のコードはlevel1でやった場合の結果です。

見ると、クラスタ1と2には外れ値がありません。どういうことかというと、クラスタ1と2は複数の下位クラスタを持っており、それらの平均をとると特定のトピックが突出しなくなってしまうのです。逆に、クラスタ3,4,5は下位クラスタがないか少ないため、あるトピックが突出して出てくる、というわけです。
同じようにlevel2の箱ひげ図を出してみましょう。

分類を細かくすると、ひとまずすべてのクラスタに外れ値が発生するようになりました。ですが、クラスタ5のように

さらに下位クラスタをもっていると、外れ値の外れ方もそれほど極端ではありません。

クラスタ2も同様です。
このように、クラスタへのトピック割り当てを箱ひげにしてやれば、どのクラスタにどのトピックが割り当てられているのか、ぱっとわかるようにできます。

各文書へのトピック割り当てを検討する

上のやり方は平均値をとっているため、そもそもそのトピックへの割り当て確率がどのような分布をしているのか、が見えません。狭い分布なのか、広い分布なのか、がはっきりしないと、平均値をとることが妥当なのか、なんとも言えないところです。そこで、同じように箱ひげ図を各トピックへの割り当て確率の分布として表現してみましょう。

# topic probability list by cluster
topicPlistL1 <- lapply(seq(1,length(levels(hierStr$level1))),function(x){
  topicProbability[][hierStr$level1==x,1:20]
})
boxplot(topicPlistL1[[1]])

として、level1のクラスタ1における各トピックへの割り当て確率の分布を箱ひげ図にしてみましょう。

横軸がトピック1～20、縦軸が各トピックへの割り当て確率です。全体的にばらつきが非常に大きくなっています。
すでにみたように、level1のクラスタ1は多数の下位クラスタをもっているため、特定のトピックへの割り当てが突出することはありませんでした。それがこの分布からも確認できます。
では、クラスター3はどうでしょうか？

トピック6への割り当て確率が高いところで分布しており、他のトピックが小さいところで分布しています。非常にはっきりと出てきており、わかりやすい形となっています。

トピックの中身は？

以上のようにして、各クラスタに固有のトピックが割り当てられていることが確認できました。では、そのように割り当てられてトピックの中身はどのようなものなのでしょうか？
LDAで抽出された各トピックを代表する単語は

termsDF <- get_terms(model,100)

で抽出することができます。このときのmodelはLDAでモデリングしたモデルです。第二引数の100は出力する単語の数なので、上記はトピック数20でLDAを行ったモデルから、各トピックごとに100単語を抽出する、という意味になります。
さて、level1のクラスタ3にはトピック6が高い確率で割り当てられているのでした。そのトピック6とは何ものかというと、

> termsDF[1:20,6]
 [1] "patient"      "chemotherapi" "group"        "rate"         "surgeri"      "outcom"       "recurr"      
 [8] "local"        "cancer"       "node"         "differ"       "p"            "resect"       "complic"     
[15] "hospit"       "lymph"        "age"          "studi"        "result"       "vs"  

化学療法（chemotherapi）、手術（surgeri）、再発（recurr）、がん（cancer）、リンパ（lymph）といった単語が出現しています。つまり、再発性のがんにたいする化学療法や手術に関する文献であることがなんとなくわかってきます。また、"p"や"studi"という単語から、治験にかかわる文献なのではないか、とも推測できます。
だいぶはっきりはしてきましたが、それでもなんとなくぼんやりした印象は変わりません。もうちょっとなんとかできないでしょうか？

代表的なセンテンスを抽出する

RにlexRankrというパッケージがあります。読み込んだ文章をセンテンスで区切り、類似度でグラフを作ったときに中心性の高いセンテンスを抽出する、というものだそうです。

www.cs.cmu.edu

そういえば、自分がテキストマイニングをやり始めのころ、同じような動作をするコードを見つけたことがありました。それをパッケージ化したもののようです。まぁ、パッケージになっているのであれば使わせてもらいましょう、ということで

install.packages("lexRankr")

でインストールし、サンプルを参考にして

library(lexRankr)
sumClust <- sapply(seq(1,5),function(x){
  testSum <- paste(dataMedDF[hierStr$level1==x,2],collapse=" ")
  testSum <- stringi::stri_split_boundaries(testSum, type = "sentence")[[1]]
  out <- lexRank(testSum,n=3)
  return(out[3][1])
})
sumClust <- unlist(sumClust)

とします。
level1のそれぞれのクラスタに属するabstractをすべて結合し、ひとつの文書とした後、それをセンテンスに分解し、中心性を計算して代表的なセンテンス上位3文ほどをとってくる、というものです。ちなみにセンテンスが増えてくると計算がかなり遅くなるので、注意してください。

さて、気になるlevel1のクラスタ3はどうだったでしょうか。

> sumClust[3]
$sentence
[1] "20.0 % (p = 1.0) for the obese patients,  respectively, without a significant difference."                                                                                                                                                                              
[2] "Significant outcome differences between groups were found in operative and anesthesia times (longer in the RS group), and in estimated blood loss, intraoperative transfusion, length of stay, and postoperative complications (all  higher in the open surgery group)."
[3] "However, there was no significant difference between the HIPEC treatment groups (groups 3, 4, 5, and 6) (p>0.05)."    

1文目と3文目は「有意な差がなかった」で2文目が「有意な差があった」と言っています。臨床試験の結果と思われますが、それぞれ反対の結論です。もちろん、試験の中身がわからないのでなんとも言えませんので、もしかしたらこのクラスタ3はもっと下の方でさらに下位カテゴリに分類できるのかもしれません。
ついでといってはなんですが、緑のbreast_cancerがいっぱいだったlevel1のクラスタ5はどうでしょうか？

> sumClust[5]
$sentence
[1] "In vitro cytotoxicity of the investigated compounds was tested on three human cancer cell lines HeLa, A549 and MDA-MB-453 using MTT assay."
[2] "pachycarpa were investigated on DU145 (prostate cancer cell line) and MCF-7 (breast cancer cell line) cell lines."
[3] "The cytotoxicity of the PSG toxin against cancer cell lines was determined using the MTT assay."

あんまり詳しくないのではっきりはわかりませんが、cytotoxicityといっているので細胞毒性に関する文献なのでしょう。ここからさらに文献タイトルを見ていけば、どんな分類になっているのかはっきりするのではないでしょうか。
※今回はそこまでやっていません。

まとめ

というわけで、LDAとその結果を用いた階層的クラスタリングによって、かなりきれいに文書の分類ができることがわかりました。また、文書への割り当て確率を箱ひげ図にすることによって主要トピックの抽出が比較的容易になること、代表センテンスの抽出でさらにそれを深掘りしていけることがわかりました。
今回はデンドログラムを指定した高さで分割していましたが、例えばデンドログラム上で任意の高さをクリックするとその点でのクラスタリング、トピック割り当ての箱ひげ図、主要ワードと主要センテンスがクラスタ毎に出力される、なんてことができたら、文献調査ももっと楽にできるんじゃなかろうかな、と考えている次第です。
こういった手法で、結果の解釈が面倒なLDAをもっと身近に扱っていけたらいいな、と思うのであります。