読者です 読者をやめる 読者になる 読者になる

CodeIQ Blog

自分の実力を知りたいITエンジニア向けの、実務スキル評価サービス「CodeIQ(コードアイキュー)」の公式ブログです。

R+igraph問題「友好関係ネットワークから派閥を検出」 @kztakemoto さんによる解説記事 #R #igraph #sna

CodeIQ中の人、millionsmileです。

大学でバイオインフォマティクスやネットワーク理論の研究をしている竹本さんから、Rを使ったソーシャルネットワーク分析問題の解説記事です!

==============================

はじめに


近年、ネットワーク分析はデータ分析の一種であり、さまざまな分野で注目されています。マーケティング施策やコンテンツ・商品推薦などの文脈で、ネットワーク分析を要求される機会はますます増えていくことだと思います。

一方で、「ネットワーク分析は取っ付きにくい」と感じる方も多いのではないでしょうか。そんなことはありません。最近では、ネットワーク分析を支援するソフトウェアやツールが充実してきています。特に、データ分析でよく用いられる"R"では、ネットワーク分析を支援するパッケージ"igraph"が利用可能です。

この問題では、R+igraphを利用してネットワーク分析で多用されるコミュニティ分析に取り組んで頂き、ネットワーク分析をより身近なものに感じて頂くことを目的としました。

ただ、このようなパッケージを用いる場合、基本的に関数の使い方に終始してしまい、分析そのものの背景がおろそかになりがちです。そこで、このエントリーでは問題解説に加え、ネットワーク分析を行う際に参考となる情報についても言及したいと思います。

参考文献はすべて閲覧可能ですので、より詳しく知りたい方はそちらも参照してください。

問題のおさらい

学校や会社などでは、いくつかの派閥ができることがしばしばです。どうやら、とある空手クラブにはふたつの派閥があるとのことです。このクラブにおける友好関係ネットワークは下図のようになっていますが、その派閥どのようになっているかは分かりません。ネットワーク分析の手法、特に、コミュニティ検出法を用いてその派閥を見つけてください。

問題の解説


Q1からQ3までの、正解率はほぼ100%でした。ただ、解答を見て気になった点がありましたので、以下に補足説明したいと思います。

ネットワークの読み込みとファイル形式

Q1. GraphML形式で記述されるZacharyの空手クラブにおける友好関係ネットワークを読み込み、ノードやエッジの属性(アトリビュート)を出力してください。 特に、ノードのラベル、エッジの重み、ノードのメンバーシップ(実際の派閥)を出力してください。

ネットワークは関数read.graphを用いて読み込むことができます。ファイル形式はGraphML(後で解説します)ですので、その指定format="graphml"を忘れないようにします。

> library(igraph) # igraphライブラリの読み込み
# graphml形式のネットワークデータの読み込み
> g <- read.graph("karate.GraphML",format="graphml")

関数summaryを用いれば、このネットワークにどのような属性が含まれるか参照することができます。

> summary(g)
IGRAPH UNW- 34 78 -- Zachary's karate club network
attr: name (g/c), Citation (g/c), Author (g/c), Faction (v/n), name
  (v/c), id (v/c), weight (e/n)

ネットワーク(グラフ)の属性は、以下のようにして参照できます。

> g$name #ネットワークのタイトル
[1] "Zachary's karate club network"
> g$Citation #ネットワークの出典元
[1] "Wayne W. Zachary. An Information Flow Model for Conflict and Fission in Small Groups. Journal of Anthropological Research Vol. 33, No. 4 452-473"

ノードとエッジの属性はそれぞれ、関数Vと関数Eを用いて参照できます。

> V(g)$name #ノードのラベル
 [1] "Mr Hi"    "Actor 2"  "Actor 3"  "Actor 4"  "Actor 5"  "Actor 6" 
 [7] "Actor 7"  "Actor 8"  "Actor 9"  "Actor 10" "Actor 11" "Actor 12"
[13] "Actor 13" "Actor 14" "Actor 15" "Actor 16" "Actor 17" "Actor 18"
[19] "Actor 19" "Actor 20" "Actor 21" "Actor 22" "Actor 23" "Actor 24"
[25] "Actor 25" "Actor 26" "Actor 27" "Actor 28" "Actor 29" "Actor 30"
[31] "Actor 31" "Actor 32" "Actor 33" "John A"  
> V(g)$Faction #ノード(人)が属する派閥(Faction)の番号
 [1] 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2
> E(g)$weight # エッジの重み(友好関係の相対的な強さ)
 [1] 4 5 3 3 3 3 2 2 2 3 1 3 2 2 2 2 6 3 4 5 1 2 2 2 3 4 5 1 3 2 2 2 3 3 3 2 3 5
[39] 3 3 3 3 3 4 2 3 3 2 3 4 1 2 1 3 1 2 3 5 4 3 5 4 2 3 2 7 4 2 4 2 2 4 2 3 3 4
[77] 4 5

GraphMLはとても便利なファイル形式です


GraphML形式での読み込みは上記のようにformat=で指定することが簡単に読み込むことができますが、GraphMLについては理解できたでしょうか。

実際のネットワークには、このように多くの属性が含まれることがしばしばです。これらの属性を含めネットワークを記述するには、どのよう形式を用いればよいのでしょうか。そのひとつとして、GraphML形式が注目されています。

GraphMLは、ネットワークに特化した拡張可能なマークアップ言語(XML)のひとつです。GraphMLを用いれば、複雑になりがちな属性付きネットワークデータをひとつのファイルにまとめて簡単に記述することが可能になります。

例えば、今回使用した空手クラブネットワークのGraphMLファイルは以下のようになっています。関数summaryを使わなくても、このファイルの中身を見れば、どのような属性が含まれているかを知ることができます。

<!-- 上略 -->
<!-- 属性の指定 -->
  <key id="name" for="graph" attr.name="name" attr.type="string"/>
  <key id="Citation" for="graph" attr.name="Citation" attr.type="string"/>
  <key id="Author" for="graph" attr.name="Author" attr.type="string"/>
  <key id="Faction" for="node" attr.name="Faction" attr.type="double"/>
  <key id="name" for="node" attr.name="name" attr.type="string"/>
  <key id="weight" for="edge" attr.name="weight" attr.type="double"/>
  <graph id="G" edgedefault="undirected">
    <data key="name">Zachary&apos;s karate club network</data>
    <data key="Citation">Wayne W. Zachary. An Information Flow Model for Conflict and Fission in Small Groups. Journal of Anthropological Research Vol. 33, No. 4 452-473</data>
    <data key="Author">Wayne W. Zachary</data>
    <!-- ノード情報 -->
    <node id="n0">
      <data key="Faction">1</data> <!-- 派閥 -->
      <data key="name">Mr Hi</data> <!-- ラベル(名前) -->
    </node>
<!-- 中略 -->
  <!-- エッジ(二項関係の)情報 -->
    <!-- n0とn1に間にエッジが描かれることを意味する -->
    <edge source="n0" target="n1">
      <data key="weight">4</data> <!-- 重み -->
    </edge>

ネットワーク描画

Q2. ネットワークを描画してください。エッジの重み(友好間関係の強さ)を描画に反映させたり、色やレイアウトに工夫を施したりしてください。

実のところ、igraphは描画に関してはあまり優秀であるとは言えません。この点に関しては、Cytoscapeの方が適しています。また、Cytoscapeもネットワーク分析に関する様々なプラグインが充実しています。ネットワーク分析にはCytoscapeや他のツールやソフトウェアの選択肢もあることを覚えておくとよいでしょう。

igraphにおけるネットワーク描画の例は以下のようです。その他のオプションもありますので、参考にしてみて下さい。

> plot(g,
+ vertex.size=15, #ノードの大きさ
+ vertex.shape="rectangle", #ノードの形
+ vertex.label=V(g)$name, #ノード属性nameをノードラベルにする。
+ vertex.color=ifelse(V(g)$Faction==1,"Pink","Lightgreen"), #ノード属性Factionを用いてノードに色づけ
+ vertex.label.color="gray50", #ノードのラベルの色
+ vertex.label.font=2, #ノードのラベルのスタイル 1: 普通, 2: 太字, 3: 斜体, 4: 太字斜体, 5: ギリシャ文字
+ vertex.frame.color="white", #ノードの枠の色
+ vertex.label.cex=0.8, #ノードラベルの文字サイズ
+ edge.width=E(g)$weight, #エッジ属性weightをエッジの太さとする
+ edge.color="gray80", #エッジの色
+ layout=layout.fruchterman.reingold) #ネットワークのレイアウト手法

ネットワーク描画の結果は以下のようになります。

コミュニティ検出

Q3. igraphの関数を用いてコミュニティ検出を行ってください。igraphには様々な手法が利用可能です。いくつかの手法を試してください(可能ならば全て試しください)。

適切な関数を用いれば正解です。ただ、何人かの解答者から「コミュニティ検出法の細かい点は良くわかっていないのですが・・・」というコメントがありましたので、以下に補足説明します。

コミュニティ検出の基礎


コミュニティ検出法の多くは以下のように定義されるQ値(Clauset et al. 2004)が最大となるようなメンバーシップを見つける問題として捉えられます。
Q=\frac{1}{2E}\sum_{ij}\left(A_{ij}-\frac{k_ik_j}{2E}\right)\delta(c_i,c_j)
ここで、A_{ij}はネットワークの隣接行列を、Eはエッジ数を、 k_iはノードiの次数(近傍ノード数)を示します。\delta(c_i,c_j)は、もしノードiとjが同じコミュニティに属している場合1を返し、そうでなければ0を返すような関数です。このQ値は「コミュニティ内のエッジ密度が偶然に得られるそれと比較してどの程度大きいか」を表す指標になっています。

ただ、ネットワークが大きい場合、全ての可能なメンバーシップ状態を探索することほぼ不可能です。そのためなんからの手法を用いてQ値が最大となるようなメンバーシップ状態を求める手法が一般に用いられます。

有名なコミュニティ検出法


igraphでは、様々な手法に基づくコミュニティ検出法が利用可能です。詳細は竹本がWebで公開する資料でご覧頂けますが、有名な手法は次のようなものが挙げられます。なお、グラフオブジェクトにエッジの重みなどが含まれている場合、それらを考慮して計算してくれます。

Girvan-Newman法(Girvan & Newman, 2002

辺媒介性に基づくカットによる分割を通してコミュニティを検出します。

> edge.betweenness.community(g)
Graph community structure calculated with the edge betweenness algorithm
Number of communities (best split): 6 
Modularity (best split): 0.345299 
Membership vector:
   Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8 
       1        1        2        1        3        3        3        1 
#下略

貪欲法に基づく手法(Clauset et al., 2004

隣接するノード/コミュニティをマージさせ、評価値の高いマージを選択することで、最適な分割を見つけます。

> fastgreedy.community(g)
Graph community structure calculated with the fast greedy algorithm
Number of communities (best split): 3 
Modularity (best split): 0.4345215 
Membership vector:
   Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8 
       2        2        2        2        3        3        3        2 
#下略

固有ベクトルに基づく手法(Newman, 2006

グラフラプラシアンを用いて最適なカットを見つけることで、コミュニティを検出します。

# ARPACK(大規模固有値問題のためのライブラリ)のオプションを与えることができます。
> leading.eigenvector.community(g,options=list(maxiter=1000000, ncv=5))
Graph community structure calculated with the leading eigenvector algorithm
Number of communities (best split): 4 
Modularity (best split): 0.3934089 
Membership vector:
   Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8 
       1        3        3        3        1        1        1        3 
#下略

焼きなまし法に基づく手法(Reichardt & Bornholdt, 2006

焼きなまし法を用いてエネルギー最小となるよう(Q値が最大になるような)なコミュニティのメンバーシップ状態を見つけることでコミュニティを検出します。

> spinglass.community(g)
Graph community structure calculated with the spinglass algorithm
Number of communities: 5 
Modularity: 0.2079796 
Membership vector:
   Mr Hi  Actor 2  Actor 3  Actor 4  Actor 5  Actor 6  Actor 7  Actor 8 
       2        2        2        2        1        1        1        2 
#下略

メンバーシップに関しては、それぞれ、$membershipで参照することができます。以下は貪欲法を用いた例です。

> comm_greedy <- fastgreedy.community(g)
> comm_greedy$membership
 [1] 2 2 2 2 3 3 3 2 1 1 3 2 2 2 1 1 3 2 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1

どの手法を選択すればよいのでしょうか?:計算時間と検出精度のトレード・オフ


さて、様々な手法が利用できるのは良いことなのですが、結局のところ、どれを使えばよいのでしょうか。それぞれの手法には利点と欠点があります。目的と状況に応じて選択することが重要です。

基準になるのは計算時間(とにかく早く計算したい)と検出精度(より最適な解をを得たい)だと思います。Fortunato (2010)はコミュニティ検出の多くの手法について、計算時間と検出精度の視点から比較しています。

計算時間では、貪欲法(fastgreedy.community)が最短です(O(n\ln^2n)、ここで nはノード数)。計算時間が短い順に挙げると次のようになります:貪欲法 < Girvan-Newman法 < 固有ベクトル法 < 焼きなまし法。

検出精度では、焼きなまし法(spinglass.community)が最高です。精度の高い順に挙げると次のようになります:焼きなまし法 > 固有ベクトル法 > 貪欲法 > Girvan-Newman法。ただ、この結果は、あるベンチマークを用いた検証から得られたものですので、解釈には注意が必要でしょう。あくまでひとつの基準として考えて下さい。

この結果から分かるように、計算時間と検出精度はトレード・オフの関係にあります。これは手法を選択するひとつの判断基準になるでしょう。

Q値によるコミュニティ検出の問題点


多くのコミュニティ検出法はQ値を評価値として用いています。しかしながら、このQ値に基づくコミュニティ抽出法には解像度限界(Fortunato & Barthelemy, 2007)という問題があることが知られています。簡単に言えばこれは、ある程度小さなコミュニティを見逃してしまう場合がある、という問題です。

この問題は、Q値が偶然に得られるコミュニティ内のエッジ密度(ヌルモデル)との比較に基づいていることが原因となっています(Fortunato & Barthelemy, 2007; Fortunato ,2010)。具体的には、小さいコミュニティが偶然に形成される確率は、大きいコミュニティのそれよりも大きくなるからです。

この問題を避けるために、Li et al. (2008)はヌルモデルを必要としない指標Modularity Density(D値)を提案しました。これは以下のように、あるモジュール内に属するノードG_iの内部次数(あるノードに注目した場合、同じコミュニティに属する近傍ノード数)の平均d_{\mathrm{in}}(G_i)と外部次数(異なるコミュティに属する近傍ノード数)の平均d_{\mathrm{out}}(G_i)の差の全て(m個)のコミュティについて合計したもの、として定義されています。
D=\sum_{i=1}^md_{\mathrm{in}}(G_i)-d_{\mathrm{out}}(G_i)
ただ、このD値を用いた場合でも、解像度限界を完全に避けることはできません。Li et al. (2008)は、更にこのD値の拡張版を提案しています。

自身でコミュニティ検出アルゴリズムを実装するさいは、評価値をQ値の代わりにD値を用いてもよいかもしれません。

コミュティ分析:検出精度の確認

Q4. コミュニティ抽出によって検出された派閥と実際の派閥(Q1で読み込んだノードのメンバーシップ)を比較し、その検出精度(一致率)を計算してください。複数のコミュニティ検出法を試して、それぞれの検出精度を算出してください。

Q4の正解率はぐっと下がります(7人中3人でした)。これはみなさん方向性は間違っていなかったのですが、以下のふたつの点について配慮されていませんでした。

問題中に記載される「派閥はふたつある」という情報が用いられていない。


コミュニティ検出法は最大Q値における分割を採用するので、コミュニティが2つになるとは限りません。例えば、貪欲法の場合、最大Q値は0.43でその時の、コミュニティ数は3つです。

> comm_greedy
Graph community structure calculated with the fast greedy algorithm
Number of communities (best split): 3 
Modularity (best split): 0.4345215 

コミュニティ数を調整するためには、マージ行列(どのようにノード/コミュニティがマージされた過程を表す行列で、$mergesで参照可能です)と任意の計算スッテプ数を入力として、関数community.to.membershipを用います。この関数は任意の計算ステップでの分割におけるメンバーシップとそれぞれのコミュニティサイズを返してくれます。なお、この関数は階層的にクラスタリングする以下のアルゴリズムのみで有効です。

  • edge.betweenness.community
  • walktrap.community*1
  • fastgreedy.community

例えば、貪欲法を例に挙げると以下のようになります。

以下のように、各ステップにおけるコミュニティの数をみると、ステップ32でコミュニティ数が2つになっていることが分かります。

> for(i in 1:nrow(comm_greedy$merges)){
+ d <- community.to.membership(g, comm_greedy$merges, steps=i)
+ cat("step ",i," #Comms=",max(d$membership)+1,"\n")}
#中略
step  31  #Comms= 3 
step  32  #Comms= 2 
step  33  #Comms= 1 

念のために、上のようにチェックしましたが、マージの過程を考えれば、コミュニティ数 = ノード数 - ステップ数、なのでステップ数は以下のようにしても求められます。

> vcount(g) - 2
[1] 32

まとめると、以下のようにして、コミュニティがふたつに分かれる場合のメンバーシップを得ることができます。

> mem_comm_greedy_c2 <- community.to.membership(g,
+ comm_greedy$merges, steps=vcount(g)-2)$membership
> mem_comm_greedy_c2
 [1] 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1

なお、焼きなまし法の場合はオプションspins=で最大コミュニティ数を指定することで、コミュニティ数を調整することができます。

> spinglass.community(g, spins=2)$membership
 [1] 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2

コミュニティ検出法から得られたメンバーシップの番号が必ずしも実際の派閥の番号対応するとは限らない。


例えば、以下のようです。

> V(g)$Faction
 [1] 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2
> mem_comm_greedy_c2
 [1] 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1

番号に基づいて一致率を計算する場合、これらの番号を揃える必要がありますが、多くの解答者はその作業を行っていませんでした。

適切に一致率を計算するためには、例えば、以下のようします。これは、全問正解者のひとりであるtonbo様のコードを参考にしています(許可を頂いて掲載しています)。

# 番号を揃える
> mem_comm_greedy_c2 <- ifelse(mem_comm_greedy_c2==0,1,2)
# 一致率の計算
> sum(ifelse(V(g)$Faction==mem_comm_greedy_c2,1,0)) / vcount(g) * 100
[1] 100

他の正解者は違った方法で一致率を計算していました。ここは様々なやり方があると思います。みなさんはどのように計算しますか?

なお、同様にして、他の手法についても一致率を計算すると以下のような結果になります。

  • edge.betweenness.community: 94.11765%
  • walktrap.community: 100%
  • spinglass.community: 100%

Fortunato (2010)が言うように、Girvan-Newman法はあまり検出精度が高いとは言えないことが分かります。

まとめ


R+igraphを用いると、いい意味で、ネットワーク分析がお手軽に行えることがお分かり頂けたと思います。このエントリーでは問題の解答に加え、ネットワーク分析について覚えておくと分析を行う際に参考となる情報についても言及しました。まとめると以下のようです。

  • ネットワーク分析を支援するツールはたくさんありますが、igraphは多くの分析手法が利用可能であることに加え、Rとの併用で柔軟な分析ができるという利点があります。
  • 多くの属性をもつネットワークを記述する場合GraphML形式が役立ちます。
  • コミュニティ検出法の多くは、なんらかの最適化手法を用いてある評価値が最大となるようなメンバーシップ(どのノードがどのコミュニティに属するか)をみつける問題として捉えられます。
  • 評価値としてよく用いられるQ値には解像度限界という問題点が存在します。いくつか別の評価値が提案されていますが、いまだに議論の余地が残っています。
  • 計算時間と検出精度にはトレード・オフの関係があります。

この空手クラブのネットワークはコミュニティ検出法の性能評価のベンチマークとしてよく使われます。もし、新たなコミュニティ検出法を提案した際には、このネットワークを使って性能評価してみてはいかがでしょうか。

また近年では、コミュニティの重複(ひとつのノードが複数のコミュニティに属すること)を考慮した手法も注目を浴びています。ここでは言及しませんでしたが、Rパッケージのひとつlinkcomm(Kalinka & Tomancak, 2011)を用いることで、分析することができます(より詳しく知りたい方はWeb資料をご覧下さい)。興味をもたれましたら、お試し下さい。

この問題の挑戦者は7名で、決して多くの方にご興味頂けたというわけではありませんでしたが、今後この解説がネットワーク分析を理解する上で、また実際に分析を行う上で、みなさまのお役に立つことを願っています。

竹本和広 @kztakemoto

==============================

いやー、竹本さんの解説記事、とても丁寧でたっぷりと贅沢に説明してくれていて素晴らしいですね。さすがホンモノの研究者!!!

記事の中に、より詳しく知りたい人向けにWeb資料がリンクされていますが、このリンク先にある「R+igraphではじめる生物ネットワーク解析」という発表資料はとてもよくできていて一読の価値アリです。竹本さんがこの資料の発表しているのを聴きにいったのですが、テンポのよいしゃべりで、発表を聴いたあと暫くは興奮が冷めやらぬ感じでした。個人的にはoverlapping communitiesの話しが面白かったです。
http://www.slideshare.net/kztakemoto/r-seminar-on-igraph

https://codeiq.jp/ace/takemoto_kazuhiro/q325
f:id:codeiq:20130626101312p:plain

今後も面白い問題をだしていきますので、CodeIQをよろしくお願いしますー!

エンジニアのための新しい転職活動!CodeIQのウチに来ない?の特集ページを見る

*1:ランダムウォークに基づくアルゴリズム([http://jgaa.info/accepted/2006/PonsLatapy2006.10.2.pdf:title=Pons & Latapy, 2006])です。