ネットワーク分析をもうちょっと勉強（実践編）
前回は 竹本先生の資料 や WEBページ を拝見し R で行うネットワーク分析について新たに知ることができた関数、分析手法についてまとめてみた。さっそく Non-overlapping 手法を実務で実践したが早くも壁にぶちあたってしまった。実務で取り扱うデータを分析する際にはつきものだ。。。今回分析したネットワークはこんな感じ。ノード数=446、エッジ数=792、図ではエッジの矢印が見えないかもしれないが有向グラフ。（前回同様サイトスケープにくわせOrganicレイアウトで表示 ※背景色のみ白に変更）
さっそくエラー多発
何も気にせずにNon-overlapping手法のコミュニティ抽出関数を実行したらエラーが多発。
> # おまじない
> ls(all=T)
character(0)
> rm(list=ls(all=TRUE))
>
> # パッケージの読み込み
> library(igraph)
>
> # ワーキングディレクトリの変更
> setwd("d:/R/igraph") # 作業ディレクトリの変更
>
> # ネットワークの読込み
> d <- read.table("エッジデータ.txt", header=T) # データフレームで読込み
> g <- graph.data.frame(d[,2:3], directed = TRUE) # データフレームからグラフオブジェクトに変換
> g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除
>
> eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割
> rw <- walktrap.community(g) # ランダムウォークに基くアルゴリズム
> fg <- fastgreedy.community(g) # 貪欲アルゴリズム Q値が高くなるいようにペアをまとめていく
以下にエラー fastgreedy.community(g) :
At fast_community.c:537 : fast greedy community detection works for undirected graphs only, Unimplemented function call
> le <- leading.eigenvector.community(g) # スペクトラル最適化法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
警告メッセージ：
In leading.eigenvector.community(g) :
At community.c:1357 :This method was developed for undirected graphs
> ml <- multilevel.community(g) # 多段階最適アルゴリズム
以下にエラー multilevel.community(g) :
At community.c:2421 : multi-level community detection works for undirected graphs only, Unimplemented function call
> sp <- spinglass.community(g) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
以下にエラー spinglass.community(g) :
At clustertool.cpp:285 : Cannot work with unconnected graph, Invalid value
> op <- optimal.community(g) # 全探索（時間がかかるため小さなネットワークのみ有効）
以下にエラー optimal.community(g) :
At glpenv05.c:70 : glp_malloc: no memory available
, CPU time exceeded
>うひょ～ なんだ？ 辺媒介性とランダムウォークは問題なし、貪欲アルゴリズムと多段階最適アルゴリズムは無向グラフのみか。（スペクトラル最適化法は警告）全探索はマシンのスペック不足のようだ。最も試したかった焼きなまし法はネットワークが分断されているということか。上の図を見ても確認できるがネットワークは既にいくつかのコミュニティに分かれているってことだ。前に igraph を使っていたときは clusters 関数を使い連結成分のみのネットワークに分割し、その単位でコミュニティ抽出していたことを思い出した。前処理、前処理。
ネットワークの連結成分を確認する
igraph の clusters 関数は連結成分を確認することができる。（Rで学ぶデータサイエンスの#8ネットワーク分析にも記載されている）
> g.cl <- clusters(g) #連結成分を確認 > g.cl$membership # 連結成分 [1] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1 5 1 [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 7 1 1 1 1 1 1 1 4 1 1 1 1 :(中略) [419] 1 1 1 1 8 8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 > g.cl$csize # 連結成分のサイズ [1] 400 9 5 12 4 2 11 3 > g.cl$no # 連結成分の数 [1] 8
そうそう、この関数で分断されているネットワークがいくつあるか、またノードはどの群に属しているか確認できる。前はこの関数を使い Excel でネットワークを分割していたが今回は R で何とかしてみるか。。。
decompose.graph 関数
分断されたネットワークを連結成分単位に分割するコードを書こうかと考えたが、ひょっとしたら 既に関数があるかもと思い約300ページある igraph のマニュアルを確認してみることにした。まずは目次から。すると decompose.graph 関数というものがあった。decompose = 分解 ということは？？？ ひょっとして発見した！？ ということでさっそく試す。（dで始まる関数だったの結構早く見つかって良かった）
> dcg <- decompose.graph(g) # 連結成分でグラフを分割する modeは"weak"みたい > length(dcg) # グラフの連結数 [1] 8 > dcg [[1]] IGRAPH DN-- 400 752 -- + attr: name (v/c) [[2]] IGRAPH DN-- 9 8 -- + attr: name (v/c) ：（中略） [[8]] IGRAPH DN-- 3 2 -- + attr: name (v/c)
ビンゴ！ ネットワークは8つに分割されていて、各ネットワークはそれぞれ配列っぽく格納されているじゃないか。これは使える！
非連結ネットワークのコミュニティを抽出する
ということで、非連結ネットワークに対してコミュニティ抽出するコードを書いてみた。私のＲ経験値ではこんなコード。（この例では焼きなまし法だが他の手法でも同様）
# おまじない
ls(all=T)
rm(list=ls(all=TRUE))
# パッケージの読み込み
library(igraph)
# ワーキングディレクトリの変更
setwd("d:/R/igraph") # 作業ディレクトリの変更
# ネットワークの読込み
d <- read.table("エッジデータ.txt", header=T) # データフレームで読込み
g <- graph.data.frame(d[,2:3], directed = TRUE) # データフレームからグラフオブジェクトに変換
g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除
# 連結成分でネットワーク
dcg <- decompose.graph(g)
# 連結ネットワーク単位にコミュニティーと中心性を計算する
sp.df.all <- c() # コミュニティ分割結果データフレーム
for (i in 1:length(dcg)){
set.seed(1) # シードを固定
sp <- spinglass.community(dcg[[i]]) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
sp.df <- cbind(i, as.data.frame(sp$names), as.data.frame(sp$membership)) # データフレーム形式に変換
sp.df.all <- rbind(sp.df.all, sp.df) # コミュニティーに結果を追加
}
colnames(sp.df.all) <- c("i", "NodeName", "Community") # 列名変更
write.table(sp.df.all, "Community.txt", row.name=F, col.names=T, sep="\t", quote=F, append=F) # ファイル出力
焼きなまし法によるコミュニティ抽出結果
上のコードを使い焼きなまし法によるコミュニティ抽出を行い冒頭のネットワーク図のノードを色分けするとこのようになった。（描画はcytoscapeにノード属性として結果ファイルを食わせて色分け）
でたぁっ 初回にしては結構いい感じ～～ ここまできたら中心性指標から特徴的なノードも見てみるか。
媒介中心性から特徴的なノードを確認
とり急ぎ媒介中心性を使い特徴となりそうなノードも確認。媒介中心性の指標値が大きいノードを赤色にしてグラデーションをかけてみた。（描画はコミュニティ同様cytoscapeにノード属性として結果ファイルを食わせて色分け）
betweenness(g, directed=TRUE, weights=NULL) # 媒介中心性
中心性指標は分断されているネットワークでも計算できる（前のバージョンでは計算できなかった関数があったような気もするが。。。）が、中心性指標によってはコミュニティ抽出同様、連結成分ごとに計算しなければいけない場合もあるかも。。。その時はコミュニティ抽出の for ループの中で計算しちゃえば良いってことだ。
中心性指標を比較してみる
分断されているネットワークの場合には指標値が異なる結果を返す中心性指標がある為、分断されていない最も大きなネットワークで比較。とりあえず相関。
i <- 1 # 分断されていない最も大きなネットワーク cl <- closeness(dcg[[i]], mode="out", weights=NULL, normalized=TRUE) # 近接中心性 dg <- degree(dcg[[i]], mode="out") # 次数中心性 ev <- evcent(dcg[[i]], directed=TRUE, weights=NULL)$vector # 固有ベクトル中心性 pr <- page.rank(dcg[[i]], directed=TRUE, damping=0.85, weights=NULL)$vector # PageRank bp <- bonpow(dcg[[i]], exponent=0.2) # ボナチッチのパワー中心性 bw <- betweenness(dcg[[i]], directed=TRUE, weights=NULL) # 媒介中心性 nc <- cbind(cl, dg, ev, pr, bp, bw) # 中心性指標まとめ library(psych) pairs.panels(as.data.frame(nc), smooth=FALSE, density=FALSE, ellipses=FALSE, scale=TRUE) # 散布図
コミュニティ抽出 焼きなまし法と辺媒介性の比較
最期にコミュニティ抽出で比較した焼きなまし法と辺媒介性の2つを比較した結果。辺媒介性ではコミュニティが138もある。（対象ノード数は400） dendPlot 関数で確認してみると、げげ～ 多くのノードがまとまった比較的大きなコミュニティーが数個と、ほぼ１つで形成されているだろうノードが数多く存在する結果となっているようであった。今回分析したネットワークの特徴？ 有向グラフだから？ ということで今回は焼きなまし法を採用した。パフォーマンスも良いということなので。
> i <- 1 # 分断されていない最も大きなネットワーク > sp <- spinglass.community(dcg[[i]]) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す > eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割 > max(sp$membership) # 焼きなまし法のコミュニティ数 [1] 16 > max(eb$membership) # 辺媒介性のコミュニティ数 [1] 138 > dendPlot(eb)
まずはNon-overlapping手法を実践。今度はOverlapping手法とか実践してみながら今のデータ分析に有効か確認しよう。
