次数分布を固定してランダムグラフを生成する
はじめに
グローバルクラスタリング係数や平均パス長、度数中心性など、ホールネットワークに関する統計量について、客観的に見てその数値がどれほど大きなものなのかを知りたい時がある。
このような場合、ヌルモデルと比較し、グラフの特徴量が統計的にどれほど逸脱しているのかを調べるという方法がある。よく使われるのは、オリジナルのグラフと同様の次数分布を持つランダムグラフと比べるというものである。
この記事では、Rのigraphパッケージを使って、次数分布を固定してランダムグラフを生成する方法について紹介する。
グラフ読み込み
今回はigraphに付属している空手クラブネットワークデータを用いる。
library(igraph) library(tidyverse) library(ggnetwork) net <- graph("Zachary")
とりあえず可視化
ggnetworkでグラフを可視化する。後ほど、ランダムグラフも可視化するが、オリジナルのグラフでのノードの座標を最初に取得しておき、使い回せるようにしておく。
hruchterman_reingold = layout.fruchterman.reingold(net) gnet <- ggnetwork(net,layout=hruchterman_reingold) g <- ggplot(gnet,aes(x=x,y=y,xend=xend,yend=yend)) g <- g + geom_edges(size=0.5) g <- g + geom_nodes(size=3,col="blue") g <- g + geom_nodelabel_repel(aes(label=vertex.names)) g <- g + ggtitle("original karate club network") g <- g + theme_blank() g
画像がこちら。なんとなく、クラスタリング係数が高く、平均パス長は大きい、ような気がする。
実際に計算してみると
> cc_original <- transitivity(net) #グローバルクラスタリング係数 > pl_original <- mean_distance(net) #平均パス長 > > cc_original [1] 0.2556818 > pl_original [1] 2.4082
さて、これらの値はどこまで特徴的なのだろうか？
ランダムグラフ生成
igraphでは、sample_degseq(in.deg) で生成できる。in.deg はそれぞれのノードの次数のベクトル。無向グラフの場合はin.degのみ、有効グラフの場合はin.degで入次数の分布、out.degで出次数の分布を指定する。
オリジナルのネットワークの次数分布を取得しておく。
degree_dist <- igraph::degree(net) g <- ggplot(data=tibble(degree=degree_dist),aes(x=degree)) g <- g + geom_histogram() g <- g + ggtitle("degree distribution of karate club network ") g
グラフの生成方法はオプション引数のmethodで指定する。igraphにはsimple (デフォルト), simple.no.multiple,vl の3種類が用意されている。
simple
ノードから適当にn本のエッジを生やして、それを結合するという方法。シンプルな方法だが、多重エッジやセルフループが発生する。つなげ合わせる2つのノードの持つエッジの相手が重複する可能性を排除していないためである。
degree_dist <- igraph::degree(net) simple_random_net <- sample_degseq(degree_dist)
こうして得られたグラフだが、simplify() して多重エッジとセルフループを取り除くと、エッジ数が少し減ることが確認できる。
> simple_random_net IGRAPH 3c731f6 U--- 34 78 -- Degree sequence random graph + attr: name (g/c), method (g/c) + edges from 3c731f6: [1] 2-- 6 1-- 4 3--14 17--28 1--34 19--28 12--33 7--18 9--29 21--33 5--34 1--14 3--33 1--28 8-- 8 22--27 25--31 [18] 7--32 9--33 1--22 19--20 11--34 1--20 2-- 7 4--34 1--30 18--34 4--32 2-- 5 10--34 1-- 1 6--14 1-- 3 9--17 [35] 6--33 8--24 5--34 16--30 4--23 1--31 33--34 3--24 25--27 33--34 2--26 1-- 3 3--33 28--34 26--33 1-- 3 23--34 [52] 4--29 14--32 10--34 2--29 9--32 11--20 34--34 3--24 11--25 4--15 3--15 1--14 2--31 1--16 7--24 2--34 3--13 [69] 21--32 6--32 24--30 1-- 2 2--34 13--33 8--30 33--33 9--34 26--31 > simplify(simple_random_net,remove.multiple = TRUE,remove.loops = TRUE) IGRAPH 40d3c01 U--- 34 65 -- Degree sequence random graph + attr: name (g/c), method (g/c) + edges from 40d3c01: [1] 1-- 2 1-- 3 1-- 4 1--14 1--16 1--20 1--22 1--28 1--30 1--31 1--34 2-- 5 2-- 6 2-- 7 2--26 2--29 2--31 [18] 2--34 3--13 3--14 3--15 3--24 3--33 4--15 4--23 4--29 4--32 4--34 5--34 6--14 6--32 6--33 7--18 7--24 [35] 7--32 8--24 8--30 9--17 9--29 9--32 9--33 9--34 10--34 11--20 11--25 11--34 12--33 13--33 14--32 16--30 17--28 [52] 18--34 19--20 19--28 21--32 21--33 22--27 23--34 24--30 25--27 25--31 26--31 26--33 28--34 33--34
simple.no.multiple
生成ロジックは simple と同じだが、こちらは多重エッジやセルフループが発生したらそのグラフは破棄し、多重グラフやループのないグラフが生成されるまで生成を続ける。したがって平均次数が多い場合は計算時間がかかり、また、可能なグラフ集合からの均一なサンプリングではない。
この方法で生成したランダムグラフを可視化する。
gnet <- ggnetwork(simple_no_multiple_random_net,layout=hruchterman_reingold) g <- ggplot(gnet,aes(x=x,y=y,xend=xend,yend=yend)) g <- g + geom_edges(size=0.5) g <- g + geom_nodes(size=3,col="blue") g <- g + geom_nodelabel_repel(aes(label=vertex.names)) g <- g + ggtitle("sampled karate club network (method: simple.no.multiple)") g <- g + theme_blank() g
オリジナルのネットワークに比べると、クラスタ化の度合いが小さく、ノード間の距離は近いように思われる。
先ほどの simple によって発生させたグラフと異なり、simplify() の前後でエッジ数は変化しない。
> ecount(simple_no_multiple_random_net) [1] 78 > ecount(igraph::simplify(simple_no_multiple_random_net,remove.multiple = TRUE,remove.loop = TRUE)) [1] 78
vl
- 難しいので説明は省略。
こちらについても不均一サンプリングになるらしい
igraphの公式ドキュメント曰く
The “vl” method is a more sophisticated generator. The algorithm and the implementation was done by Fabien Viger and Matthieu Latapy. This generator always generates undirected, connected simple graphs, it is an error to pass the in.deg argument to it. The algorithm relies on first creating an initial (possibly unconnected) simple undirected graph with the given degree sequence (if this is possible at all). Then some rewiring is done to make the graph connected. Finally a Monte-Carlo algorithm is used to randomize the graph. The “vl” samples from the undirected, connected simple graphs unformly. See http://www-rp.lip6.fr/~latapy/FV/generation.html for details.
有意性判定
実際に先程のクラスタリング係数と平均パス長をヌルモデルからの乖離によって評価してみる。vlメソッドによってランダムグラフを1000個発生させ、そこから計算した指標の分布の中でオリジナルネットワークの指標がどのあたりに位置するのかを検討する。
sample_cc <- c() sample_pl <- c() for(i in 1:1000){ vl_random_net <- sample_degseq(degree_dist,method = "vl") sample_cc <- c(sample_cc,transitivity(vl_random_net)) sample_pl <- c(sample_pl,mean_distance(vl_random_net)) }
各指標の分布は以下のようになる。
#cc hist g <- ggplot(data = tibble(cc = sample_cc),aes(x = cc)) g <- g + geom_histogram(binwidth = .01) g <- g + geom_vline(xintercept = cc_original,col = "red") g #pl hist g <- ggplot(data = tibble(pl = sample_pl),aes(x = pl)) g <- g + geom_histogram(binwidth = .01) g <- g + geom_vline(xintercept = pl_original,col = "red") g
これらの分布が正規分布に従っていると仮定（平均パス長に関してはずいぶん右に歪んでいるが…）して、平均と標準偏差とz値を計算し、p値を算出する。
cc_mean = mean(sample_cc) pl_mean = mean(sample_pl) cc_sd = sd(sample_cc) pl_sd = sd(sample_pl) cc_zscore <- (cc_original - cc_mean) / cc_sd pl_zscore <- (pl_original - pl_mean) / pl_sd
> 1-pnorm(cc_zscore) [1] 0.1035129 > 1-pnorm(pl_zscore) [1] 0.0003767265
ランダムグラフから得られた指標の分布と比較すると、オリジナルのネットワークは両指標とも統計的に逸脱した値であることが分かる。
さいごに
ここで紹介した一連の方法は、あるグラフの中でどのような構造が特徴的なのかを探索する際に使う程度に留めておいたほうがよいかもしれない。つまり、他のグラフと比較する際に上の統計的な逸脱度を比べる、ということはしないほうが良いかもしれない。
というのも、クラスタリング係数や平均パス長の異なるグラフ間の比較方法について検討した論文Comparing Brain Networks of Different Size and Connectivity Density Using Graph Theoryによれば、これらの指標はノードサイズNや平均次数kの影響を大きく受けるため、元に異なるサイズや密度のグラフを生の数値によって比較することは困難である（こうしたバイアスはN-independentやk-independent等と呼ばれている）。
今回紹介したようなヌルモデルベースのメトリックのほうがむしろNやkの影響を受けてしまうようで、安易な使用には注意が必要である。
