はじめに

ネットワーク分析に興味を持ち、共立出版のRで学ぶデータサイエンスシリーズの「ネットワーク分析」を読みました。
適用範囲は未知数ですが、なかなか面白いと思いました。
２,3回に分けて内容を簡単にまとめたいと思います。
初めはネットワークの構造やノードの特徴を表す指標についてです。

igraphとggraph、tidygraphパッケージを中心に使用していきます。
こちらで掲載しているコードの詳細はgithubにあげています。

ネットワーク分析について

ネットワークとは

ネットワークとは頂点とそれらをつなぐ辺で構成された頂点との関係を表現するものです。
例えば、頂点は人、辺を人々の繋がりであるとするとネットワークは組織内のコミュニケーション関係を表現するものになります。
インターネットも、頂点はWebページ、辺はリンクであり、一種のネットワークであると言えます。
また、グラフ理論においては、頂点をノード、辺をエッジとも呼ばれるため、この記事では、そちらを用いたと思います。

[出典]

何ができるのか

SNS上の人々のコミュニケーションの繋がりを見て、コミュニケーションの中心的な人物を抽出したり、企業間の関係を見てどのような企業が協業しているのかを抽出、路線図から最短経路を見つけ出すようなことができます。

[出典]

ネットワークデータの加工

ネットワークデータの加工は少し他のデータと異なります。
ネットワークデータ(GLMフォーマット)そのものを読み込めば、比較的簡単に加工ができます。
今回は１から入力してグラフデータを作成しようと思います。
初めにグラフのノードとエッジの関係を表す隣接行列を作成します。
$i$番目のノードから$j$番目のノードへのエッジが存在することを、$i$行$j$列に1を立てることにより表現します。
0はそのノード間にエッジが存在していないことを意味しています。
なお、エッジの方向性があるネットワークを有向グラフ、方向性を考えない(エッジの全てが双方向に存在)場合のネットワークを無向グラフと呼びます。

library(igraph)
library(tidygraph)

A <- matrix(c(
  0,1,1,1,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0, 
  1,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,0,1,0,0, 
  1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0, 
  1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0, 
  0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0, 
  1,0,0,0,1,0,1,0,0,0,1,0,0,0,1,0,1,0,0,0, 
  0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0, 
  0,1,0,0,0,0,0,0,1,1,0,1,0,0,0,0,0,0,1,1, 
  0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,1, 
  0,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,0,0,
  0,1,1,1,0,1,0,0,0,0,0,1,0,0,0,1,0,0,0,0, 
  1,0,1,0,0,0,0,1,0,0,1,0,1,1,0,0,0,1,0,0, 
  1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0, 
  0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0, 
  0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0, 
  0,0,0,0,1,0,1,0,0,0,1,0,0,1,0,0,1,0,0,0, 
  0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0, 
  0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,1, 
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1, 
  0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0),
  nrow = 20, ncol = 20, byrow = TRUE)

colnames(A) <- c(1:20)
rownames(A) <- c(1:20)

matrixからグラフデータに変換します。
内部では、どのノードからどのノードにエッジが存在するのかを示す情報が記載されています。

> g1 <- graph_from_adjacency_matrix(A)
> g1
IGRAPH 3f6c467 DN-- 20 86 -- 
+ attr: name (v/c)
+ edges from 3f6c467 (vertex names):
 [1] 1 ->2  1 ->3  1 ->4  1 ->6  1 ->12 1 ->13 2 ->1  2 ->8  2 ->11 2 ->13 2 ->18 3 ->1  3 ->11
[14] 3 ->12 4 ->1  4 ->11 5 ->6  5 ->7  5 ->16 5 ->17 6 ->1  6 ->5  6 ->7  6 ->11 6 ->15 6 ->17
[27] 7 ->5  7 ->6  7 ->15 7 ->16 8 ->2  8 ->9  8 ->10 8 ->12 8 ->19 8 ->20 9 ->8  9 ->10 9 ->18
[40] 9 ->20 10->8  10->9  10->14 11->2  11->3  11->4  11->6  11->12 11->16 12->1  12->3  12->8 
[53] 12->11 12->13 12->14 12->18 13->1  13->2  13->12 14->10 14->12 14->16 15->6  15->7  15->17
[66] 16->5  16->7  16->11 16->14 16->17 17->5  17->6  17->15 17->16 18->2  18->9  18->12 18->19
[79] 18->20 19->8  19->18 19->20 20->8  20->9  20->18 20->19

tidygraphパッケージで提供されている、グラフデータをtidyに扱うtbl_graph型にも変換します。
tbl_graph型で扱うことによりdplyrとの相性がよく処理を記述することができます。

> g1_tbl <- as_tbl_graph(g1, direted = F)
> g1_tbl
# A tbl_graph: 20 nodes and 86 edges
#
# A directed simple graph with 1 component
#
# Node Data: 20 x 1 (active)
  name 
  <chr>
1 1    
2 2    
3 3    
4 4    
5 5    
6 6    
# … with 14 more rows
#
# Edge Data: 86 x 2
   from    to
  <int> <int>
1     1     2
2     1     3
3     1     4
# … with 83 more rows

ネットワーク描写

ネットワークの構造や関係性を見ていくためには数値で見るよりも、可視化した方が格段にわかりやすいです。
plot関数でも描写はできますが、おしゃれではないです。
ggplot2のラッパーであるggraphを利用するとおしゃれなグラフを描けることができます。

library(ggraph)

# グラフ描写
g1_tbl %>%  
  ggraph(layout = "kk") +
  geom_edge_link(alpha=0.8, colour = "lightgray") + 
  scale_edge_width(range = c(0.1,1)) +
  geom_node_point(aes( size = 4)) +
  geom_node_label(aes(label = name),repel = TRUE)

ネットワークの指標

それではネットワークの指標についてです

密度

密度(density)はネットワークにおいて貼ることができる全ての辺に対して、実際に存在する辺の数の割合です。
頂点の数を$n$、グラフにある辺の数を$m$とすると、密度は次のように表すことがきます。

density = \frac{2m}{n(n-1)}

そして、　Rでは次のように求めることができます。

> graph.density(g1_tbl)
[1] 0.2263158

推移性

ある頂点$i,j,k$の三点が相互に繋がっている場合、この関係を推移的と言います。
例えば、友達の友達は自分の友達でもあるような関係です。
推移性(transitivity)は、ネットワークにおいて推移的な三点がどれだけあるのかを表す指標です。
Rでは次のように求めます。

> transitivity(g1_tbl)
[1] 0.28125

相互性

有向辺は3種類あり、双方向の辺と一方向のみの辺(2種類)です。
相互性(reciprodcity)は、有効辺において双方向の辺がどれだけあるのかを表す指標です。
例えば、あるグループにおいて相互にコミュニケーションにある人が多ければ高くなり、コミュニケーションが活発でも一方向のコミュニケーションしかなければ、低くなります。

数式で表すと次のようになります。

reciprodcity = \frac{a}{a+c+d}

ここで、$a$は双方向の辺の辺の数、$c,d$は一方向の辺の数を表しています。
Rでは簡単に求めることがきます、今回のグラフは初めから全ての辺が有向辺を仮定していたので、１となります。

> reciprocity(g1_tbl)
[1] 1

ネットワーク中心性

中心性とは、ネットワークにおける各頂点の重要性を評価したり、比較するための指標です。
tidygraphではここに記載するだけではなく、様々な中心性指標抽出のコマンドが提供されています

次数中心性

次数とは、各ノードに接続しているエッジの数のことです。
次数中心性は、この次数つまり、エッジが多く存在しているノードにおいて大きくなる指標です。

# 描画関数
plot_tbl_data <- function(tbl_data, title = "centrality"){
  tbl_data %>%  
    ggraph(layout = "kk") +
    geom_edge_link(alpha=0.8, colour = "lightgray") + 
    scale_edge_width(range = c(0.1,1)) +
    geom_node_point(aes( size = centrality)) +
    geom_node_label(aes(label = name),repel = TRUE)+
    ggtitle(title)
}

# 次数中心性
g1_tbl.deg <- g1_tbl %>% 
  mutate(centrality = centrality_degree())

plot_tbl_data(g1_tbl.deg, "centrality:degree")

隣接中心性

隣接中心性(closeness centrality)とは、他の頂点との距離が小さい頂点ほど、より中心的であるとする指標です。
隣接中心性は次のように計算することができる。

C(i) = \frac{1}{\sum_{j=1}^n d_{ij}}

$d_{ij}$は頂点$i$から$j$への最短距離の合計である。

# 隣接中心性 
g1_tbl.cls <- g1_tbl %>% 
  mutate(centrality = centrality_closeness())

plot_tbl_data(g1_tbl.cls,"centrality:closeness")

媒介中心性

対象となるノードを削除した時、全体のネットワークが２つに分離するようなノードをブリッジと呼びます。
社会ネットワークにおいて、その人がいなくなると組織が二分されるような人がこの媒介中心性が高くなる人です。
媒介中心性は、このブリッジを高く評価する指標です。

# 媒介中心性

g1_tbl.bet <- g1_tbl %>% 
  mutate(centrality = centrality_betweenness())

plot_tbl_data(g1_tbl.bet, "centrality:betweenness")

個人的には、次数中心性が高いものが媒介中心性が高い気がします。

固有ベクトル中心性

単純に多くのノード接続していることのよりも、接続しているノードの中心性が高い場合の方がそうでない場合よりも重要性が高いと言えそうです。
固有ベクトル中心性は、無向グラフの隣接行列における第一固有ベクトルを指標としたものです。
ある無向グラフの隣接行列を$A=(a_{ij})$とし、そこにお希る頂点の中心性を成分とする列ベクトルを$c=(c_j)$とすると、頂点iの中心性$c_i$は次のように表すことができる。

c_i = \frac{1}{\lambda}　\sum_{j=1}^n a_{ij} c_j

ここで$\lambda$は比例定数で正の値とする。

# 固有ベクトル中心性
g1_tbl.eig <- g1_tbl %>% 
  mutate(centrality = centrality_eigen())

plot_tbl_data(g1_tbl.eig, "centrality:eigen")

Page Rank

固有ベクトル中心性では、無向グラフのみ扱えましたが、有向グラフにまで適用できる中心性指標です。

# page rankの算出
g1_tbl.pr <- g1_tbl %>% 
  mutate(centrality = centrality_pagerank())

plot_tbl_data(g1_tbl.pr, "centrality:page rank")

情報中心性

情報中心性(infomation centrality)は、ネットワークに含まれる頂点間の最短経路以外の経路や、経路の長さも考慮した中心性指標です。
点間の各経路の長さの逆数を重みづけし、その総和を規格化することで求める。

# 情報中心性
g1_tbl.inf <- g1_tbl %>% 
  mutate(centrality = centrality_information())

plot_tbl_data(g1_tbl.inf, "centrality:infomation")

コミュニティ抽出

ネットワークのノードのクラスタリング=コミュニティ抽出を行います。
igraphでもコミュニティ抽出の処理が行えるが、今回はtidygraphを用います。
tidygraphではここに記載するだけではなく、様々なコミュニティ抽出のコマンドが提供されています。

媒介中心性によるコミュニティ抽出

ネットワーク内で相対的に密度の高いサブグループを求めていきます。
ネットワーク内での密度の高さを凝縮性(cohesion)と呼び、次式で求めることができる。

cohesion=
\frac{
\frac{\sum_{i \in S}\sum_{j \in S}a_{ij}}{n_s(n_s-1)}}{
\frac{\sum_{i \in S}\sum_{j \notin S}a_{ij}}{n_s(n-n_s)}
}

ここで、$n$はネットワーク全体の頂点数、$n_s$はサブグループ$S$に属する頂点数である。
この式の分子はサブグループないの頂点同士の関係の密度であり、分母はサブグループないの頂点とサブグループ会の頂点との関係密度を表現している。
この値が１より大きい場合、サブグループは外部と区別可能な凝縮性を持っていると言えます。
このサブグループの探索には、エッジに対して算出した媒介性中心性指標を用います。
中心性が高いエッジの削除を繰り返すことによりグラフを分割し、凝縮性の高いコミュニティを抽出します。

# 描画関数
plot_tbl_data <- function(tbl_data, title){
  tbl_data %>%
    ggraph(layout = "kk") +
    geom_edge_link(alpha=0.8, colour = "lightgray") + 
    scale_edge_width(range = c(0.1,1)) +
    geom_node_point(aes(colour = nclass, size = 4)) +
    geom_node_label(aes(label = name), repel = TRUE) +
    ggtitle(title)
}

# 媒介中心性によるクラスタリング
g1_tbl.bet_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_edge_betweenness()))
plot_tbl_data(g1_tbl.bet_C,"betweenness")

固有ベクトルに基づくコミュニティ抽出

# 固有ベクトル中心性によるクラスタリング
g1_tbl.eig_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_leading_eigen()))
plot_tbl_data(g1_tbl.eig_C,"eigenvector")

ランダムウォークに基づくコミュニティ抽出

# ランダムウォークによるクラスタリング
g1_tbl.walk_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_walktrap()))
plot_tbl_data(g1_tbl.walk_C,"random walk")

その他のコミュニティ抽出

上記に示したコミュニティ抽出法の他にも様々な方法があります。
簡単に紹介だけします。

# 情報中心性によるクラスタリング
g1_tbl.inf_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_infomap()))
plot_tbl_data(g1_tbl.inf_C,"infomap")

# スピングラス法によるクラスタリング
g1_tbl.spin_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_spinglass()))
plot_tbl_data(g1_tbl.spin_C,"spin glass")

# グリーディアルゴリズムによるクラスタリング
g1_tbl.gr_C <- g1_tbl %>% 
  mutate(nclass = as.factor(group_fast_greedy()))
plot_tbl_data(g1_tbl.gr_C)

ネットワークの類似度

次はネットワークの類似度についてです。
ここまでのネットワークと比べるために次のようなネットワークBを作成しました。

ハミルトン距離

ハミルトン距離は、2つのネットワークの異なるエッジの数になります。

# ハミルトン距離
> # ハミルトン距離
> sum(abs(A-B))
[1] 24

相関係数

二つのネットワークにおける相関を算出します。
ネットワークの相関は２つ考え方があります。
一つ目は２つのネットワークにおいて隣接行列の同じ場所にある成分同士を比較する方法です。
これはノードの入れ替えを考えないものです。
もう一つは２つのネットワークにおいて隣接行列のノードを相関が高くなるように、入れ替えものです。
詳しい、相関係数の算出方法はこちらをご覧ください。
Rではぞれぞれ次のように表すことができます。

> # ネットワーク同士の相関(入れ替えなし)
> gcor(A,B)
[1] 0.8212116
> # ネットワーク同士の相関(入れ替えを考える)
> gscor(A,B)
[1] 0.8212116

中心化の類似性

複数のネットワークで、どのノードが中心的かを求め比較することができます。
特にネットワークが時間と共に変化していく場合、ネットワークの変化を観察することができるようになります。
この分析を中心化共鳴分析と呼び、多変量解析のコレスポンデンス分析を隣接行列に適用する分析です。
MASSパッケージに含まれているコレスポンデンス分析の関数correspを利用します。

３つのグラフを比較したいため、今までの二つに加えて新しいネットワークCを作成しました。

# 中心性の類似性 中心化共鳴分析
library(sna)
library(MASS)

degree_df <- data.frame(A_deg = degree(A),B_deg = degree(B),C_deg = degree(C))
cre <- corresp(degree_df, nf = min(dim(degree_df)) - 1)

# データの形成とグラフ化
cre$rscore <- cre$rscore %>% as.data.frame %>% mutate(type = "nord")
cre$cscore <- cre$cscore %>% as.data.frame %>% mutate(type = "graph")
cre_df <- rbind(cre$rscore, cre$cscore) %>% 
  mutate(name = c(1:20,"A","B","C"))

p <- ggplot(cre_df,aes(x = V1, y = V2, label = name)) +
  geom_hline(yintercept = 0, colour = "gray75") +
  geom_vline(xintercept = 0, colour = "gray75") +
  geom_text(aes(colour = type), alpha = 0.8, size = 5) 
p

Aでは17、Bでは10、cでは19,16が中心的であると言えそうです。

各軸の累積寄与率を求めてみます。
二軸で100%となっているため、上記のグラフにおいてデータを全ての情報を表現できていると言えます。

> eig <- cre$cor^2
> cumsum(100 * eig / sum(eig))
[1]  59.38814 100.00000

コア/周辺構造の抽出

ネットワークないで相互的に密に結合したノードが存在しており、それらをネットワークのコアと呼ぶます。
また、コア以外のノードを周辺と呼びます。
複数のネットワークにおいて、コアを比べることによりネットワークの重要な違いや/変化を知ることができます。

コアとなるノードは相互に完全に結合している場合が理想出来である。
しかし、一般的には完全に結合している場合は少なく、そこまでをコアとすればいいのかわわからない。
そこでいくつかのノードをコアと仮定しそれらが相互に完全に結合したネットワークを作成する。
このネットワークと対象のネットワークとの相関係数が最大となるコアの組み合わせを探索することで、コアの抽出を行う。
この探索を遺伝的アルゴリズムにより行います。
遺伝的アルゴリズムのパッケージであるGAを用います。
今回はAとCのコアの違いを観察します。

library(GA)

make.pat.matrix <- function(x) {
  pat.matrix <- matrix(NA, length(x), length(x))
  pat.matrix[which(x == 1), which(x == 1)] <- 1
  pat.matrix[which(x == 0), which(x == 0)] <- 0
  diag(pat.matrix) <- NA
  return(pat.matrix)
}

# Aのコア抽出
matrix.cor.A <- function(x) {return(gcor(A, make.pat.matrix(x)))}
ga.core.A <- ga(type = "binary", fitness = matrix.cor.A,nBits = nrow(A))
summary(ga.core.A)

A_core = t(ga.core.A@solution) %>% as.data.frame()
g1_tbl <- g1_tbl %>% 
  mutate(nclass = as.factor(A_core$V1))
plot_tbl_com_data(g1_tbl,"estimete core by GA")

# Cのコア抽出
matrix.cor.C <- function(x) {return(gcor(C, make.pat.matrix(x)))}
ga.core.C <- ga(type = "binary", fitness = matrix.cor.C,nBits = nrow(C))
summary(ga.core.C)

C_core = t(ga.core.C@solution) %>% as.data.frame()
g3_tbl <- g3_tbl %>% 
  mutate(nclass = as.factor(C_core$V1))
plot_tbl_com_data(g3_tbl,"estimete core by GA")

ネットワークAのコア

ネットワークCのコア

大きく違いますね。
次数が多いからといってコアになるとは限らないようです。

おわり

以上で、ネットワークの指標編を終わります。
次は統計的ネットワーク分析についてをまとめようと思います。