【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で計算してグラフ化するとこんな感じになります。

f:id:wanko_sato:20170723113002p:plain

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")

f:id:wanko_sato:20170723122801p:plain

・・・なんのこっちゃわかりません。
傾向があるのかどうかも不明だし、かといって全部のグラフを一個一個みていくのはちょっと、あ、いや、かなり嫌です。
ちなみに出現頻度の低い単語を落としても、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でクラスタ数を指定し、所定の箇所でデンドログラムを切断します。その結果得られた各クラスターに属するデータの平均値を算出してグラフにします。一枚のグラフに重ねて表示するとわけわからないことになってしまったので、画面分割での表示です。

f:id:wanko_sato:20170723124451p:plain

はい、出ました。
クラスタ数を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")
})

f:id:wanko_sato:20170723124752p:plain

ぐしゃっといろんなのが混じっていたり、数が多すぎてごちゃごちゃになっていたりするクラスタがあります。クラスタ数9では、まだ十分に分離できていないようです。
ということで、numHを増やしていったグラフを出してみましょう。

numH=15の場合。
f:id:wanko_sato:20170723124930p:plain
f:id:wanko_sato:20170723124940p:plain
そこそこきれいに分離できていそうだし、特徴もしっかり見えてきそうです。が、なんとなくもうちょいな感じ。

numH=20の場合。
f:id:wanko_sato:20170723125117p:plain
f:id:wanko_sato:20170723125129p:plain
これくらいだとちゃんと分離できて、かつ特徴もしっかり捉えられていそうな感じ。
何度もいうようですが、結果の解釈はこの分野に詳しい方と協議すべき、と思います。

単語の出力

グラフだけ出しても、じゃあ各クラスタにはどんな単語が入っているんだ?と突っ込まれること間違いなしなので、

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といった単語が見えます。突発性疾患、造血系の疾患への応用を試みた文脈でしょうか。とはいえ、今回はがんを検索ワードに含めていたため、これらの単語は頻出しないのが正しいあり方と感じます。

というような感じで

単語の出現パターンをクラスタリングして特徴を抽出する、というのは個人的になかなか面白い手法と考えています。後はここから得られた結果をどう解釈していくか、というところをしっかり詰めていけば、もっと良い分析になりそうな予感がします。

おまけ

クラスタ分割の可視化

f:id:wanko_sato:20170723131234p:plain

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)

とすると、こんな絵が描けます。

f:id:wanko_sato:20170723131614p:plain

あ、いや、デンドログラムをigraphで書き直しただけなんですけどね。

個々のノードを追いかけていくと、どの段階で分割されたか、がほんのりわかるようになっています。
例えば一番右端あたりは3回分割された後、ずっとそのまま分割されないで維持されています。hclustでやっていることを考えると、もともとしっかり分離されたクラスタである、つまり特徴がはっきりしていて他と区別しやすい特徴をもったデータである、ということが言えるわけです。
そういう指標も併せて考えると、クラスタリングの結果ももっとわかりやすくなってくるかな、なんて思ったりします。

t-SNEをやってみる

blog.albert2005.co.jp

最近知ったもので、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)

f:id:wanko_sato:20170723132716p:plain

3次元プロットするとこんな感じになりまして、距離の離れている塊はちゃんと別のクラスタになるようになっています。
t-SNE自体はクラスタリングをしていないので、各色はnumH=20でクラスタリングした場合のものなんですが、ちゃんとそのクラスタリングにあった形での配置になっています。
赤いところ、つまりクラスタ1に相当する箇所ですが、この結果を見ると確かに分離しにくい形をしていて、なるほど確かに、という感じがします。

なにか機会があったら、t-SNEについても取り上げてみたいと思っております。

まとめ

というわけで、単語の出現パターンをクラスタリングする、という試みはなかなかうまくいきました。クラスタリングの手法をどうするか、とか、もとになる距離行列には何を使うか、といった細かいテクニカルな話は残っているにせよ、十分実用に耐えうるものじゃないかな、と考えております。後は、手動でクラスタ数を変更していたところを、たとえばShinyでインタラクティブに動作させられるようにする、とか使い勝手の部分でしょうか。

そういう話はまたいずれ。