【Rで自然言語処理】単語の出現パターンをクラスタリングする。
前回、前々回とトピックモデルに関する話題を扱いました。
wanko-sato.hatenablog.com
wanko-sato.hatenablog.com
トピックモデルは嫌いじゃない、むしろ好きなんですが、結果の解釈が面倒なのと、やっぱり確率分布を使うため、結果にどこかふわっとした感じが残るあたり、どうにもモヤってしまうのです。もちろん、そういうのを織り込み済みであることを分かってもらえれば良いのですが、LDAはどうにも説明がややこしいところがございまして、簡便な方法で、細かい説明をしなくてもすぐにわかってもらえる指標ってないのかなぁ、とあれこれ考えているところでございます。
その中でも、最近思いついた中で個人的にヒットしているのが「単語の出現パターン」でございまして。今日はそのお話をしようと思います。
そもそもの経緯
論文のabstractから研究トレンドを見つけることができないか?というのがそもそもの動機です。
研究トレンドを対象にするわけですから、経時データを扱うことになります。
トピックモデルを使って、トピックの経時推移を、とも思ったのですが、結果の解釈が面倒である上、経時推移を扱うとなるとさらにややこしくなってしまいそうです。
であるならば、もっと単純な指標、例えば単語の出現頻度の傾向から何か知見が得られたら良いのじゃないか。
早速やってみよう
対象とするデータ
今回も同じようにPubMedのデータを使います。
ただし、なるたけ経時推移がはっきりしているものが良いので、iPS細胞に関する論文を対象とすることにしました。
また、単純にiPS細胞だけにすると論文数が膨大になってしまうため、"ips stem cell cancer"で検索することとしました。
※iPS幹細胞を癌化させてがん細胞の研究に使う、というような話題がここ最近(?)のいつだったかあったため。
検索でヒットしたのは491件、うち、今回の対象となるのはabstractの情報が含まれている473件です。
PubMedの検索結果をMEDLINE形式でダウンロードし、データフレームにして分析する、という流れは前々回と同じです。MEDLINE形式のデータをデータフレームに変換するRコードは前々回の記事にありますので、そちらを参考になさってください。
データの概要
dataMedList <- lapply(dataMed,function(x){ out <- data.frame(title=x$TI, abst=ifelse(is.null(x$AB),"",x$AB), publicationDate=x$DP, year=as.numeric(strsplit(as.character(x$DP)," ")[[1]][1])) return(out) }) dataMedDF <- do.call(rbind,dataMedList) dataMedDF <- dataMedDF[dataMedDF$abst!="",] yearPnum <- table(dataMedDF$year) plot(yearPnum)
MEDLINE形式を読み込ませたdataMedから必要なデータをlapplyで抽出してリスト形式にします。今回は出版年の情報も必要なため、DP(=Publication Date)も取得します。ただし、DPは年月または年で示されているため、年の情報だけ抜き出す必要があります。
その後、do.callでデータフレームにし、abstractがブランクの文献は除を外します。
で、各年毎の論文数をtableで計算してグラフ化するとこんな感じになります。
2007年まではあまり論文がでていませんが、2008年から急激に件数が増えてきています。
2012年まで増え続けた後、件数が落ち着いてきているようです。
iPS幹細胞に関するニュースを調べれば、いつ、どんなことが起こったかがわかるのでしょうが、それは置いといて、出版の状況としてはそういうことなのだなぁ、とほんのり頭に置いといて、次に進むことにします。
単語の出現頻度を計算する
次に、各abstractを前処理してDTM(Document Term Matrix)を作成し、TDF(Term Document Frequency)を計算します。
ちなみに、使用するパッケージは""tm","text2vec","textmineR"です。インストールしていない方はどうぞ。
library(textmineR) library(text2vec) library(tm) allTexts <- as.character(dataMedDF$abst) # create key data abstID <- as.character(seq(1,nrow(dataMedDF))) abstTitle <- as.character(dataMedDF$title) abstCategory <- (dataMedDF$year) abstLength <- sapply(allTexts,function(x){ length(strsplit(x," ")[[1]]) }) # 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")
ここまでが前処理。句読点の削除、stop wordsの削除、数字の削除、ステミングを行っています。
※今回は考えませんでしたが、場合によっては数字が重要な単語に含まれていることもあるため(例えばコドンの何番、とか)、数字の削除はやらない方が良いこともあるでしょう。
# 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) # word vectorize vectorizer <- vocab_vectorizer(vocab) #create DTM dtm <- create_dtm(it, vectorizer) # Term frecency TDF <- TermDocFreq(dtm)
DTMを作るところは定型処理です。ngramをやるときはcreate_vocablaryのところでオプション指定します。で、作ったDTMからTDFを算出します。
> head(TDF[order(TDF$term_freq,decreasing=T),],20) term term_freq doc_freq idf cell cell 3712 463 0.02136833 ip ip 1264 367 0.25373354 stem stem 1117 425 0.10700622 pluripot pluripot 782 390 0.19294865 reprogram reprogram 564 201 0.85579048 express express 562 248 0.64566664 induc induc 515 368 0.25101245 differenti differenti 501 238 0.68682471 human human 483 218 0.77460033 cancer cancer 472 176 0.98861139 use use 469 271 0.55697657 gene gene 438 191 0.90682196 generat generat 394 229 0.72537338 factor factor 346 201 0.85579048 tumor tumor 304 122 1.35507434 studi studi 297 211 0.80723726 patient patient 280 114 1.42289694 develop develop 271 178 0.97731184 potenti potenti 263 189 0.91734837 embryon embryon 254 186 0.93334871
頻出上位4語はidf(Inverse Document Frequency)の値が小さいので全文書集合を通してよく出現する単語であることがわかります。一方で、その下のreprogramやhumanはidfの値が相対的に大きいので特定の文書にしか出現しない、特徴語である、といえます。
が、ここからわかるのはあくまで全体的な特徴であり、今回目標としている経時的な情報はここからは得られません。そういうわけで、ひと手間かけてあげる必要があるのです。
単語の出現パターン分析
出現頻度を経時データにする
上で計算したTDFは文書集合全体のものでした。そこで、ここから、TDFを経時データに直していきます。
今回は各論文の出版年もデータとして保持したので、それを合わせてデータを作っていくことにします。
# term frequency by year yearVec <- names(yearPnum) tdfYear <- lapply(seq(1,length(yearVec)),function(x){ i <- as.numeric(yearVec[x]) cat(i,"\n") if(nrow(dataMedDF[dataMedDF$year==i,])==1){ out <- data.frame(term=colnames(dtm), term_freq=dtm[dataMedDF$year==i,], doc_freq=ifelse(dtm[dataMedDF$year==i,]==0,0,1)) }else{ out <- TermDocFreq(dtm[dataMedDF$year==i,]) } return(out) })
はい、これでOKです。
textmineRパッケージに入っている関数TermDocFreqを、各年でフィルタしたdtmに適用していくループを回しているだけです。
このままですとリスト形式で扱いづらいので、データフレームに直します。
# create termFreq by year termFlist <- lapply(seq(1,length(tdfYear)),function(x){ if(x==1){ tdfYear[[x]][,1:2] }else{ tdfYear[[x]][,2] } }) termFDF <- do.call(cbind,termFlist) colnames(termFDF) <- c("term",yearVec)
こちらが単語出現頻度の経時データのデータフレームになります。
が、今回は使いません。
というのも、Term Frequencyは同じ文書に複数、その語が出現していてもカウントするためです。ひとつの文書にたくさんその語が出現していれば頻度が高くなってしまい、語の出現パターンとしてはちょっと扱いづらいものになってしまいます。
そこで、今回は「その語が出現した文書の数」を指標として使うことにします。
# create docFreq by year docFlist <- lapply(seq(1,length(tdfYear)),function(x){ if(x==1){ tdfYear[[x]][,c(1,3)] }else{ tdfYear[[x]][,3] } }) docFDF <- do.call(cbind,docFlist) docFDF <- data.frame(docFDF,sum=rowSums(docFDF[,2:ncol(docFDF)])) colnames(docFDF) <- c("term",yearVec)
やってることは同じですね。ただ、ここからもうひと手間かけます。
というのも、すでに見たように、各年毎の出版論文数は数件~70件と大きく異なっています。そうなると、単語の出現件数はその年の論文数に引っ張られてしまうことになるため、指標としてはいまいち、ということになってしまいます。そこで、「その語が出現した論文数」を「その年の論文数」で割った値、つまり「出現割合」を指標として採用することにします。
# calc doc fraction by year yearPnumF <- c(yearPnum,sum(yearPnum)) docFrac <- lapply(seq(1,nrow(docFDF)),function(x){ out <- docFDF[x,2:ncol(docFDF)] / yearPnumF return(data.frame(out)) }) docFracDF <- do.call(rbind,docFrac)
はい、これで出来上がり、と言いたいところですが、ひとつ問題があります。
> table(TDF$doc_freq) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 3069 762 375 254 151 134 103 79 67 58 54 34 44 23 35 33 36 32 20 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 25 20 13 23 16 13 12 9 17 13 11 7 9 7 10 6 7 7 5 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 4 7 5 5 4 5 2 7 2 4 2 9 4 5 1 6 3 1 3 58 59 60 61 62 64 65 66 67 68 69 70 71 72 73 76 77 78 79 6 3 2 2 5 1 1 5 1 3 3 5 2 3 3 2 2 1 4 80 82 84 85 86 87 88 89 91 92 93 94 95 96 97 98 100 102 103 1 1 1 1 1 1 2 1 5 1 1 1 2 1 2 3 1 1 2 104 105 106 107 111 114 116 117 118 122 124 125 128 131 132 134 135 145 150 3 2 2 1 1 2 1 1 1 1 1 1 2 1 3 2 1 1 1 164 170 176 178 179 186 189 191 201 211 218 229 238 248 271 367 368 390 425 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 463 1
今回対象としている473件の論文中、1回しか出現しない語が3069語、2回しか出現しない語が762語、となっています。全体の出現単語の種類が5806種ですから、1回2回程度しか出現しない単語が全体のトレンドに影響を与えるとは考えにくいです。ということで、今回はひとまず
docFracDFP <- docFracDF[docFDF[,18]>nrow(dataMedDF)*0.05,]
として、5%より多い論文に出現している単語のみを対象とすることにしました。
※この5%はあくまで任意です。これによって本当は重要な単語を見逃しているかもしれませんし、逆にノイズを拾っているかもしれません。そのあたりは検索した分野に詳しい方と相談して決めるのが良いでしょう。
どんなデータなのか見てみる
さて、ひとまずすでに算出した各単語の出現パターンをグラフにしてみましょう。
matplot(t(docFracDFPP[,1:16]),type = "l")
・・・なんのこっちゃわかりません。
傾向があるのかどうかも不明だし、かといって全部のグラフを一個一個みていくのはちょっと、あ、いや、かなり嫌です。
ちなみに出現頻度の低い単語を落としても、360語が残っています。それを一個一個眺めるのは意味のない苦行と思います。
そこでクラスタリングの出番。
※なんでもかんでもクラスタリングするのが良い、とは思いません。が、データ分析の初期段階で、データの性質を大づかみにとらえるのには非常に有用な手法と思っています。
早速クラスタリングしよう
# calc distance matrix edist <- dist(docFracDFP[,1:16]) # execute hclust hclustRes <- hclust(edist,method="ward.D2")
距離行列を算出し、ウォード法で階層的クラスタリングをします。
対象が360個あって、それを階層的クラスタリングするのはちょっとどうなんだ?と思わなくもないですが、明示的にクラスタ数を決めなければならないk-means法に比べるとお手軽だし、どの段階でどういう風にクラスタが別れていくのかがわかるので、階層的クラスタリングをよく使っています。
以上、クラスタリングまで終了したので、どんな具合になったかを可視化してみます。
# set cluster number numH <- 9 # cut dendrogram hclustS <- cutree(hclustRes,numH) docFracDFPP <-cbind(docFracDFP,hclustS) # create timeSeries data set timeSeries <- lapply(seq(1,length(unique(hclustS))),function(x){ rowMeans(t(docFracDFPP[docFracDFPP[,18]==x,1:16])) }) timeSeries <- do.call(rbind,timeSeries) # plot multi-graph layout(matrix(1:numH,ncol=3)) sapply(seq(1,numH),function(x){ matplot(x=yearVec,y=timeSeries[x,],ylim=c(0,1),type="l") })
numHでクラスタ数を指定し、所定の箇所でデンドログラムを切断します。その結果得られた各クラスターに属するデータの平均値を算出してグラフにします。一枚のグラフに重ねて表示するとわけわからないことになってしまったので、画面分割での表示です。
はい、出ました。
クラスタ数を9にした場合、お分かりの通り、右下のように一貫して出現頻度の高い単語、右中のように2005年あたりから増え始める単語、右上のように最初はいっぱい出てたけど最近あまり見かけない単語、などのようになんとなく分けていくことができます。ちなみに、各クラスタに属するデータを平均せずにそのまま出すとこんな感じ。
# plot multi-graph for each term layout(matrix(1:numH,ncol=3)) sapply(seq(1,numH),function(x){ matplot(t(docFracDFPP[docFracDFPP[,18]==x,1:16]),type = "l") })
ぐしゃっといろんなのが混じっていたり、数が多すぎてごちゃごちゃになっていたりするクラスタがあります。クラスタ数9では、まだ十分に分離できていないようです。
ということで、numHを増やしていったグラフを出してみましょう。
numH=15の場合。
そこそこきれいに分離できていそうだし、特徴もしっかり見えてきそうです。が、なんとなくもうちょいな感じ。
numH=20の場合。
これくらいだとちゃんと分離できて、かつ特徴もしっかり捉えられていそうな感じ。
何度もいうようですが、結果の解釈はこの分野に詳しい方と協議すべき、と思います。
単語の出力
グラフだけ出しても、じゃあ各クラスタにはどんな単語が入っているんだ?と突っ込まれること間違いなしなので、
numMember <- sapply(seq(1,numH),function(x){ length(rownames(docFracDFPP[docFracDFPP[,18]==x,])) }) termMember <- lapply(seq(1,numH),function(x){ rownames(docFracDFPP[docFracDFPP[,18]==x,]) })
として、各クラスタのメンバー数と各クラスタに属する単語を確認します。
> numMember [1] 130 10 15 24 42 27 15 3 19 6 13 13 6 5 12 5 6 3 2 4
メンバー数はこのようになっており、クラスタ1のメンバーがやたら多いことに気づきます。それ以外は割とばらけてそう。
改めて出現パターンのグラフを眺めると、クラスタ1(左上)は各データにこれといった特徴がないのか、分けづらい特徴なのか、データがぐしゃぐしゃ集まった状態になってしまっています。ならばクラスタ数を多くすればいいか、というと逆に分ける必要のないものまで無理やり分離させてしまうような感じがして、それはちょっとどうなの?と思ったりもします。
話がそれました。
次、各クラスタの単語です。
> termMember[[18]] [1] "ip" "stem" "cell"
たとえば、numH=20のときのクラスタ18、上のグラフでいうと右端の上から2番目のグラフになります。このグラフは今回取得した論文全体を通して一貫して出現している語を表しています。そのクラスタに属する単語が上の3語です。そもそも"ips stem cell cancer"で検索しているので、これらの単語が一貫して出現するのは当然です。
> termMember[[17]] [1] "hematopoiet" "transplant" "syndrom" "signific" "idiopath" "clinic"
グラフ右上、2005年あたりから急に出現割合が減った単語です。idiopathic,hematopoieticといった単語が見えます。突発性疾患、造血系の疾患への応用を試みた文脈でしょうか。とはいえ、今回はがんを検索ワードに含めていたため、これらの単語は頻出しないのが正しいあり方と感じます。
というような感じで
単語の出現パターンをクラスタリングして特徴を抽出する、というのは個人的になかなか面白い手法と考えています。後はここから得られた結果をどう解釈していくか、というところをしっかり詰めていけば、もっと良い分析になりそうな予感がします。
おまけ
クラスタ分割の可視化
hclustした後にデンドログラムをplotするのはごくごく自然な流れと思います。が、今回のようにデータがたくさんあると、デンドログラムもわかりにくいものになってしまいますし、ぱっとみてどのクラスタがどういうふうに分割されていっているのか、どうにもわかりづらいなぁ、と思っておりました。
そこで、もうちょいこれをどうにかできないか、とちょっと試してみました。
cutH <- c(rev(hclustRes$height)[1:100],0) hierStr <- cutree(hclustRes,h=cutH) hierStr <- as.data.frame(hierStr) hierStrUniq <- unique(hierStr[,1:100]) hierUniqList <- lapply(seq(1,100),function(x){ sapply(seq(1,nrow(hierStrUniq)),function(i){ paste(x,":",hierStrUniq[i,x]) }) }) hierUniqDF <- lapply(seq(2,nrow(hierStrUniq)),function(x){ cbind(hierUniqList[[x-1]],hierUniqList[[x]]) }) hierPairDF <- unique(do.call(rbind,hierUniqDF)) library(igraph) g <- graph.data.frame(hierPairDF,directed = F) lay <- layout.reingold.tilford(g,root=1) V(g)$size <- 2 V(g)$label.cex <- 0.1 layout(matrix(1:1,ncol=1)) plot(g,layout=lay)
とすると、こんな絵が描けます。
あ、いや、デンドログラムをigraphで書き直しただけなんですけどね。
個々のノードを追いかけていくと、どの段階で分割されたか、がほんのりわかるようになっています。
例えば一番右端あたりは3回分割された後、ずっとそのまま分割されないで維持されています。hclustでやっていることを考えると、もともとしっかり分離されたクラスタである、つまり特徴がはっきりしていて他と区別しやすい特徴をもったデータである、ということが言えるわけです。
そういう指標も併せて考えると、クラスタリングの結果ももっとわかりやすくなってくるかな、なんて思ったりします。
t-SNEをやってみる
最近知ったもので、t-SNEを今回のものに適用したらどうなるのかな?という、本筋とは離れた、おまけの話です。
ちなみに、t-SNEについては上記リンク先がわかりやすいので参考にしてください。
で、今回のデータでt-SNEをやるとどうなるか。
※"tsne"パッケージをあらかじめインストールしておいてください。
colVec <- rainbow(max(as.numeric(as.character(hierStr$level20))))[as.numeric(as.character(hierStr$level20))] tsneRes <- tsne::tsne(edist,k=3) rgl::plot3d(tsneRes,col=colVec,size=6)
3次元プロットするとこんな感じになりまして、距離の離れている塊はちゃんと別のクラスタになるようになっています。
t-SNE自体はクラスタリングをしていないので、各色はnumH=20でクラスタリングした場合のものなんですが、ちゃんとそのクラスタリングにあった形での配置になっています。
赤いところ、つまりクラスタ1に相当する箇所ですが、この結果を見ると確かに分離しにくい形をしていて、なるほど確かに、という感じがします。
なにか機会があったら、t-SNEについても取り上げてみたいと思っております。