ネットワーク分析をもうちょっと勉強

　マーケティング施策や商品開発の糸口を発見する為に購買履歴データでネットワーク分析を行っているのだが、今年に入りノードやエッジの数が多い隣接行列を取り扱うようになってきた。といってもまだまだビッグデータとはいいがたいが。。。cytoscapeにくわせOrganicレイアウトで表示するとこんな感じです。（他の設定はデフォルト）

今後はもっとデータ量が多くなりそうなので、ネットワーク分析をもうちょっと勉強。

　これまで参考にしていたのはRで学ぶデータサイエンスの#8ネットワーク分析。R+igraphで中心性指標を計算したり、コミュニティの抽出などはさくさくできる。会社の先輩に話をすると、もっと勉強になる資料があるよと九州工業大学の竹本先生の資料を教えてくれた。2013/3/8に開催された第2回 Rでつなぐ次世代オミックス情報統合解析研究会で発表された資料らしい。タイトルは「R+igraphではじめる生物ネットワーク解析」です。拝見させて頂くとRで学ぶデータサイエンスシリーズには記載されていないigraphの使い方、linkcommというパッケージを知ることができた。個人的に新しく知ることができた内容はこんな感じ。

・データフレームからigraphで扱うグラフに変換できる（逆もできる）
・くわせたデータの多重エッジや自己ループを削除できる
・最短パスを見つけ描画できる
・コミュニティ抽出にはスペクトル最適化法や焼きなまし法など他にも手法がある
・ネットワークの比較（和、積、差）を行うことができる
・コミュニティ抽出にはオーバーラッピングコミュニティーがある（linkcommパッケージ）

これはさっそく試してみるしかない！

個人的に新たに知ったigraphパッケージの関数

　まずはデータフレーム形式で読込みから。データは竹本先生のページにあった演習用データをダウンロードさせて頂きました。URLはこちらです。R seminar on igraph - supplementary information - Kazuhiro Takemoto

library(igraph)
d <- read.table("eco_EM+TCA.txt") # データフレームで読込み
g <- graph.data.frame(d, directed = FALSE) # データフレームからグラフに変換
g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除

df <- get.data.frame(g) # グラフオブジェクトからデータフレームに変換

感動！　これでRDBから抽出したデータをすぐ食わせられる。RODBC使っちゃおう。また、その逆もできそうだ。cytoscapeにくわせるエッジリストも簡単じゃん。ちなみに今までだとこんな感じでRにくわせていた。（RDBからデータを抽出しExcelとか使いながら整形）

g <- graph(c(1,2,1,3,1,4,1,5,1,9,2,3,2,4,3,4,5,6,5,7,5,9,6,7,6,8,7,8),n = 9,directed = FALSE)

次は最短のパスを見つけ描画してみよう。最短パスを知る関数はRで学ぶデータサイエンスにもあったが描画の工夫はなかったと思う。
例は2とのノード"D-Glucose"、"2-Oxoglutarate"間の最短パスみたいだ。

sum <- get.all.shortest.paths(g, "D-Glucose", "2-Oxoglutarate", mode="out")
V(g)[sum$res[[1]]]
E(g)$color <- "grey"
E(g, path=sum$res[[1]])$color <- "red" # ※最短経路が複数ある場合は1つ目
E(g, path=sum$res[[1]])$width <- 3
plot(g, vertex.label=V(g)$name, vertex.size=12, layout=layout.fruchterman.reingold)

いい感じ。（この例ではsum$resには2つの結果があったので1つめの結果を描画）

次はコミュニティ抽出法。はじめの2つは知っていたけど他にもこんなにありました。

eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割
rw <- walktrap.community(g) # ランダムウォークに基くアルゴリズム
fg <- fastgreedy.community(g) # 貪欲アルゴリズム Q値が高くなるいようにペアをまとめていく
le <- leading.eigenvector.community(g) # スペクトラル最適化法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
ml <- multilevel.community(g) # 多段階最適アルゴリズム
sp <- spinglass.community(g) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
op <- optimal.community(g) # 全探索（時間がかかるため小さなネットワークのみ有効）

$membershipの最大数（コミュニティ数）を比較してみるとこんな感じ。

eb	rw	fg	le	ml	sp	op
4	5	5	5	4	6	6

手法が違うのでコミュニティ数も変わるのは当然だが、資料によると計算時間とパフォーマンスはトレードオフになっているとのこと。ふむふむ。焼きなまし法（アニーリング）なんかいいのかな？

焼きなまし法で抽出したコミュニティ毎に色を変え描画するとこんな感じになった。大きなネットワークで行ったらどうなるのだろう。美しさならcytoscapeなどのソフトに食わせるしかないのかな？それは今度。

V(g)$color <- sp$membership
plot(g, vertex.size=12, vertex.label=V(g)$name, layout=layout.fruchterman.reingold)

コミュニティの状況はデンドログラムでも確認できるのだが、dendPlot関数を使うこともできるようだ。dendPlotを使うと綺麗な罫線で囲んでくれるから見やすい。fastgreedy.communityの結果はこんな感じ

dend <- as.dendrogram(fg)
plot(dend) # 罫線は出ません

dendPlot(fg) # 罫線が出るデンドログラム

※焼きなまし法だとNot a fully hierarchical community structureといったエラーになる。完全な階層になっていないの? 階層構造になっていないので表示できない。

ネットワークの比較

　竹本先生のページR seminar on igraph - supplementary information - Kazuhiro Takemotoを最期まで読むと、最期にネットワーク比較なるものが。。。igraphでできたの～って感じです。
g1、g2　といったグラフに対して行う場合にはこのように使います。

graph.union(g1, g2) # 和
graph.intersection(g1, g2) # 積
graph.difference(g1, g2) # 差 ※ただし |g1| > |g2|

igraphのマニュアルには、graph-operators に記載されていました。目次には union とかで記載されていません。300ページ近くもあるからただでさえ把握困難なのに。。。
ノード名で行いたい場合は～～.by.name関数がありました。

graph.union.by.name(g1, g2) # 和
graph.intersection.by.name(g1, g2) # 積
graph.difference.by.name(g1, g2) # 差

差の結果はこんな感じです。データは適当に作っております。
plot(g1)

plot(g2)

plot(graph.difference.by.name(g1, g2))

こんな図では確認しにくいですが g1 - g2 ができていました。g1からg2にあるエッジが消える感じです。隣接行列の演算をしているのでしょう。これは時系列変化などでの活躍間違いなしだ。

オーバーラッピングによるコミュニティー抽出

　最後にコミュニティーはオーバーラップしているとい考え方を取り入れた分析手法。linkcomパッケージを使うのか。これは初めてのパッケージ。4種類ある手法の2種類が実装されているみたいだ。さっそくlinkcomパッケージをインストールして確認。

まずは、LinkCommという階層クラスタリング手法

library(linkcomm)
g <- read.table("eco_EM+TCA.txt")
lc <- getLinkCommunities(g)

デンドログラムがでた。

階層構造と分割密度なるものが確認でき、lc$pdmax を見ると最大は0.833のようだ。
任意の密度でネットワーク図を描画するにはこんな感じ

lc <- newLinkCommsAt(lc, cutat=0.84) # cutstに任意の分割を指定
lc$clusters # クラスタリング結果を確認
lc$nodeclusters # ノード名でクラスタリング結果を確認
plot(lc, type="graph") # ネットワーク図を描画

階層クラスタリング手法の結果はこんな感じでcytoscape用エッジ属性をファイルに出力できるようだ。cytoscapeも勉強不足。

linkcomm2cytoscape(lc, interaction="pp", ea="edge_atrib.ea") # cytoscape用エッジ属性の出力

オーバーラッピングで使える関数は１つOCG(Overlapping Cluster Generator)という手法。これは貧欲アルゴリズムでLinkCommよりも良いパフォーマンスを示すとのこと。

ocg <- getOCG.clusters(g)
plot(ocg, type="graph") # ネットワーク図を描画

でたぁっ。階層クラスタリング手法でも同様だがノードの色がすごい。いくつのコミュニティーに属するかでノードが円グラフみたいになっている。これはcytoscapeでもできるのかな。。。今はそんなことはいいか。後で考えよう。

また、複数コミュニティーに属するノードはテーブルで描画することもできる

plot(ocg, type="members") # 複数コミュニティーに属するノードを見る

すげ～～～～～～～～～。竹本先生どうもありがとうございました。m(__)m　本当に感謝です！
先生の資料を教えてくれた先輩もありがとう！
自分なり今までよりネットワーク分析ができそうな気になってきた。バックグラウンドや関数の特徴などを理解しながらさっそく実務でいろいろ使ってみよう。CRANにあるigraphのマニュアルって300ページもあるから後回しにしていたのも反省。やっぱりマニュアル読まないとダメだ。

参考
・Rで学ぶデータサイエンスの#8ネットワーク分析

・R+igraphではじめる生物ネットワーク解析

・竹本先生のページ
R seminar on igraph - supplementary information - Kazuhiro Takemoto