相関行列からネットワークへ

　(昨日の記事から色々追加・修正したので、新しく書き直しました。)
　FisherのZ変換を施された相関行列でネットワーク解析を行う。

お品書き
1. FisherのZ変換の逆作業をして生の相関行列を作る
2. 有意な相関係数のみ残す
3. 相関行列からsemi-metricな距離行列を作成する
4. semi-metricな距離行列から重み付き隣接行列を作成する

1. FisherのZ変換の逆作業をして生の相関行列を作る
　某解析ソフトは相関行列を計算した結果にFisherのZ変換を施すために値が[-1, 1]から溢れる。邪魔くさいので元の相関係数に戻したい。
　入力のcsvファイルを逆関数にぶちこむだけのシンプルなもの。

matrix_Z2R.py

--------------------------------------------------

#Usage

#python3 matrix_Z2R.py [input filename] [output filename]

#

#Input

#[input filename]: csv file (correlation coefficient matrix converted Z-score by Fisher's method)

#[output filename]: output filename

#

#Output

#[output filename]: raw correlation coefficient matrix

import sys

import numpy as np

filename=sys.argv[1]

output=sys.argv[2]

matrix=np.loadtxt(filename, delimiter=',')

matrix=(np.exp(2*matrix)-1)/(np.exp(2*matrix)+1)

np.savetxt(output, matrix, delimiter=',')

--------------------------------------------------

2. 有意な相関係数のみ残す

　プログラムどうこうというより、無相関検定の復習になった。

Rmat_significnat.py

--------------------------------------------------

#Usage

#python3 Rmat_significant.py [input filename] [output filename] [number of sample] [significance level]

#

#Input

#[input filename]: Input csv filename (correlation coefficient matrix)

#[output filename]: Output csv filename

#[number of sample]: Number of samples to calculate each correlation coefficient

#[significance level]: significance level (usually 0.05)

#

#Output

#[output filename]: adjacency matrix filled with 0 or 1

import sys

import numpy as np

import scipy.stats

filename=sys.argv[1]

output_filename=sys.argv[2]

n_sample=np.int(sys.argv[3])

alpha=np.float(sys.argv[4])

#load inputfile

matrix=np.loadtxt(filename, delimiter=',')

matrix[np.isnan(matrix)]=0

#get T-value matrix by correlation coefficient matrix

Tmat=np.abs(matrix)*pow(n_sample-2, 1/2)/((1-matrix**2)**(1/2))

#get T threshold for significance

T_sig=-1*scipy.stats.t.ppf(q=[alpha/2], df=n_sample-2)

#set unsignificant values 0

matrix[Tmat<T_sig]=0

#set minus R 0

matrix[matrix<0]=0

#Output

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

3. 相関行列からsemi-metricな距離行列を作成する

　話が前後するのだが、4で距離を重みにするときに単に「距離の逆数を重みとする」のではなく、「距離に1を足した数の逆数を重みとする」のが個人的にすごいグッときたんすよね。以前遊んでいた時に、ゼロの逆数が無限大になるのがとても扱いにくかったので、なるほどなぁと。

　そんなわけで距離行列作成では「重み(相関係数)の逆数から1を引く」作業をしている。

tnorm_r2len.py

--------------------------------------------------

#Usage

#python3 tnorm_r2len.py [input filename] [output filename]

#Input

#

#[input filename]: input weighted matrix csv filename

#[output filename]: output filename

#

#Output

#

#[output filename]: t-norm length matrix

import sys

import numpy as np

inf=np.inf

input_filename=sys.argv[1]

output_filename=sys.argv[2]

matrix=np.loadtxt(input_filename, delimiter=',')

matrix[matrix==0]=inf

matrix=matrix**(-1)

matrix-=1

matrix[matrix<0]=inf

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

4. semi-metricな距離行列から重み付き隣接行列を作成する

　現時点では有意とならない相関となったノード間には結合が存在しないが、結合が存在しないノード間にも重みを見出したいとする。そこで、距離行列からダイクストラ法でノード間の最短距離を求め、Dombi t-normで変換したものを重みとしているものがあって、なるほどなぁと思いそれを採用。なお対角成分はゼロにセットしなおすので、2分グラフは未対応。

tnorm_len2wei.py

--------------------------------------------------

#Usage

#python3 tnorm_len2wei.py [input filename] [output filename]

#

#Input

#[input filename]: length matrix csv file

#[output filename]: save weighted matrix as the argv

import sys

import numpy as np

import networktools

input_filename=sys.argv[1]

output_filename=sys.argv[2]

matrix=np.loadtxt(input_filename, delimiter=',')

matrix=networktools.shortest_path_length(matrix)

matrix=(matrix+1)**(-1)

for i in range(matrix.shape[0]):

 matrix[i, i]=0

np.savetxt(output_filename, matrix, delimiter=',')

--------------------------------------------------

　そんなわけで上4つを順に呼べば、某ソフトで出力される相関行列からネットワークの距離行列と重み付き隣接行列を作成することができるようになった。

PS
　久々に文字サイズがそろわない。大量のHTMLが突っ込まれているのはわかるんだが、どの作業の時に突っ込まれるのか全くわからん…

スモールワールドネス雑考

(1) ネットワークのリンク率の違いは、Crand, Lrandの分布に影響があるか。

　ノード数は200とし、リンク率を[0.5, 0.6, 0.7, 0.8, 0.9]と変えた時、それぞれ500回ずつ作成したERグラフのclustering coefficient, characteristic path lengthをヒストグラムとして表示した。

　画像に入力し忘れたのだが、上からリンク率小さい順。

　はっきりとわかるような違いはなさそうである。

(2)ノード数が増えた時、安定した分布を得るにはERグラフをより多く作る必要があるか。

　ノード数は[30, 60, 100, 200]、ERグラフ生成回数100回、リンク率0.5の時のそれぞれのヒストグラム。

　そういう目でみれば、ノード数が増えるにしたがって若干てっぺんが尖がらなくなっているようにも見えるがはっきりしない。
　同じようにノード数[30, 60, 100, 200]で、ERグラフ生成回数30回にして再度実行。

　なんとも言えないなぁ。
　一方で、ノード数を200に固定して、ERグラフ生成回数を[30, 100, 200, 300, 500]と変えていくとやはり尖がっていく。(300～500間は変化わからんけど)

　なので再現性高めるにはそれなりにたくさんのERグラフを作成した方がよさそうだけど、実際にはそこまで影響なさそうだなという印象。

(3) CrandもLrandもいらないのでは？
　結果を眺めていて気付いたこととして、ノード数が大きく異なっても、分布の頂点が来る値ってあまり変わってなかったんですよね。例えばノード数[30, 60, 100, 200]でERグラフ生成回数500回、リンク率0.5の時の結果が次の通り。

　Crandはどれも0.5程度、Lrandはどれも1.5程度のあたりに頂点が来る。他のリンク率の場合も見てみる。

リンク率0.4

リンク率0.6

リンク率0.7

リンク率0.8

リンク率0.9

　リンク率をpとすると、Crandは大体p付近で、Lrandは2-pくらいになっている。当然ながらリンク率が1の場合にはCrand=1, Lrand=1だし、リンク率0の場合にはCrand=0になるわけなので、感覚的にはそんな変な話ではない。ただしLrandについては、リンク率が減る時に線形に減っていくわけがないので、相転移から遠いとかそういった何らかの条件付きでの話ですな。(そもそも今回のコードでは、グラフが分裂した場合は破棄している)

　前回と今回の記事と合わせて、ランダムネットワークは「どんなパラメタでもある程度同じような値を出してくる」「でもばらつきはあるよね」という雰囲気がある。はたしてこれをスモールワールドネス計算時に用いる必要があるのかとても疑問である。実質乱数みたいなもんじゃん。代わりにスモールワールドネスはランダムネットワークとか使わずCnet/Lnetとシンプルなものにしとこうや。解釈は後付けで(えー

スモールワールドネス計算時にERグラフを何個作るか

　スモールワールドネス(以後SW)は次式で定義される。

　　SW=(Cnet/Crand)/(Lnet/Lrand)
　　　Cnet: 扱っているネットワークのクラスター係数
　　　Crand: 扱っているネットワークと同じパラメタを持つERグラフのクラスター係数
　　　Lnet: 扱っているネットワークの平均最短距離長
　　　Lrand: 扱っているネットワークと同じパラメタを持つERグラフの平均最短経路長
　　　ERグラフ：Erdos-Renyi modelにより作られるネットワーク(ノード数nとリンク率p)

　Crand, Lrandが解析的に求められれば良いのだが不明。全通りの平均を求めれば確実だが、ノード数が多いと現実的な時間では終わらない。そのためErdos-Renyi modelでランダムネットワークを複数個作成してCとLを求め、その平均を用いることになる。
　では一体何個の平均をとればよいのか。
　まず1例をみてみる。個人的にこの先の未来で扱うデータはノード数10～100個程度ではないかと思うので、ひとまずノード数は80個としてみる。リンク率はとりあえずの0.5。比較するERグラフは100個でいいや。(適当)

-----------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
inf=np.inf

def shortest_path_length(lengthmat):
    n=lengthmat.shape[0]
    shortest_path_length=np.zeros((n, n))
    for i_start_node in range(n):
        i_length=np.zeros(n)+inf
        candidate=np.zeros(n)+inf
        candidate[i_start_node]=0
       
        while np.min(candidate)!=inf:
            i_length[np.argmin(candidate)]=np.min(candidate)
            candidate=i_length.reshape(n, 1)+lengthmat
            candidate=np.min(candidate, axis=0)
            candidate[i_length!=inf]=inf
       
        shortest_path_length[i_start_node]=i_length
       
    return shortest_path_length

def characteristic_path_length(lengthmat, glo=True):
    n=lengthmat.shape[0]
    shortest_path_length_mat=shortest_path_length(lengthmat)
    characteristic_path_length=np.sum(shortest_path_length_mat, axis=0)/(n-1)
   
    if glo==True:
        characteristic_path_length=np.mean(characteristic_path_length)
       
    return characteristic_path_length

def clustering_coefficient(weightedmat, glo=True):
    triangles=triangle(weightedmat)
    k=np.sum(wei2unwei(weightedmat), axis=0)
    k*=(k-1)
    k[k==0]=1
    clustering_coefficient=2*triangles/k
   
    if glo==True:
        return np.mean(clustering_coefficient)
   
    return clustering_coefficient

def wei2unwei(mat):
    n=mat.shape[0]
    adjacency=np.zeros((n, n))
    adjacency[mat>0]=1
    return adjacency   
   
def triangle(weightedmat):
    n=weightedmat.shape[0]
    triangles=np.zeros(n)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                temp=weightedmat[i, j]*weightedmat[j, k]*weightedmat[k, i]
                triangles[i]+=pow(temp, 1/3)
 
    triangles/=2
    return triangles

def ERgraph(N, p):
    graph=np.zeros((N, N))
    i=0
    while i<N-1:
        j=i+1
        while j<N:
            if np.random.uniform(0, 1)<p:
                graph[i, j]=1
               
            j+=1
        i+=1
       
    graph=graph+graph.T
    return graph

n_nodes=80
p=0.5
cc_ergraph=[]
cpl_ergraph=[]

for i in range(100):
    temp=ERgraph(n_nodes, p)
    cc_ergraph.append(clustering_coefficient(temp))
    temp[temp==0]=inf
    cpl_ergraph.append(characteristic_path_length(temp))
   
plt.hist(cc_ergraph)
plt.xlabel('clustering coefficient')
plt.ylabel('frequency')
plt.title('cc_ergraph')
plt.show()
plt.hist(cpl_ergraph)
plt.xlabel('characteristic path length')
plt.ylabel('frequency')
plt.title('cpl_ergraph')
plt.show()
-----------------------------------------------

　ちょいと足りないくらいかな、という印象。
　次に、ノード数とER生成回数を変えてどうなるかを調べてみる。ノード数は50, 100個とし、それぞれに対してERグラフの生成回数を10, 30, 100, 200, 300, 500としたときの、Cと Lの値をヒストグラムで出してみる。また正規性の確認としてシャピロウィルク検定を行う。

-----------------------------------------------
N_NODES=[50, 100]
p=0.5
n_exams=np.array([10, 30, 100, 200, 300, 500])

for n_nodes in N_NODES:
    for n_exam in n_exams:
        cc_ergraph=[]
        cpl_ergraph=[]

        for i in range(n_exam):
            cpl_candidate=inf

            while cpl_candidate==inf:
                temp=ERgraph(n_nodes, p)
                cc_candidate=clustering_coefficient(temp)
                temp[temp==0]=inf
                cpl_candidate=characteristic_path_length(temp)

            cc_ergraph.append(cc_candidate)
            cpl_ergraph.append(cpl_candidate)

        fig=plt.figure(figsize=(16, 8))
        ax1=plt.subplot(121)
        ax2=plt.subplot(122)
       
        ax1.hist(cc_ergraph)
        ax1.set_xlabel('clustering coefficient')
        ax1.set_ylabel('frequency')
        ax1.set_title('cc_ergraph ('+str(n_nodes)+' nodes, '+str(n_exam)+' exams)')

        ax2.hist(cpl_ergraph)
        ax2.set_xlabel('characteristic path length')
        ax2.set_ylabel('frequency')
        ax2.set_title('cpl_ergraph ('+str(n_nodes)+' nodes, '+str(n_exam)+' exams)')
       
        plt.savefig(str(n_nodes)+'nodes_'+str(n_exam)+'exams.png')
       
        cc_temp=scipy.stats.shapiro(cc_ergraph)
        cpl_temp=scipy.stats.shapiro(cpl_ergraph)
        print('Shapiro-Wilk test (n_nodes='+str(n_nodes)+', n_exam='+str(n_exam)+')\ncc_ergraph: p='+str(cc_temp[1])+'\ncpl_ergraph: p='+str(cpl_temp[1])+'\n')
    
-----------------------------------------------
Shapiro-Wilk test (n_nodes=50, n_exam=10)
cc_ergraph: p=0.37558773159980774
cpl_ergraph: p=0.4100293517112732

Shapiro-Wilk test (n_nodes=50, n_exam=30)
cc_ergraph: p=0.14942386746406555
cpl_ergraph: p=0.8486078977584839

Shapiro-Wilk test (n_nodes=50, n_exam=100)
cc_ergraph: p=0.7921257615089417
cpl_ergraph: p=0.3289673328399658

Shapiro-Wilk test (n_nodes=50, n_exam=200)
cc_ergraph: p=0.039140429347753525
cpl_ergraph: p=0.4688810706138611

Shapiro-Wilk test (n_nodes=50, n_exam=300)
cc_ergraph: p=0.6537032723426819
cpl_ergraph: p=0.025576237589120865

Shapiro-Wilk test (n_nodes=50, n_exam=500)
cc_ergraph: p=0.22154365479946136
cpl_ergraph: p=0.8461586236953735

Shapiro-Wilk test (n_nodes=100, n_exam=10)
cc_ergraph: p=0.7797001600265503
cpl_ergraph: p=0.38541314005851746

Shapiro-Wilk test (n_nodes=100, n_exam=30)
cc_ergraph: p=0.6236131191253662
cpl_ergraph: p=0.6676943898200989

Shapiro-Wilk test (n_nodes=100, n_exam=100)
cc_ergraph: p=0.22231294214725494
cpl_ergraph: p=0.08100102841854095

Shapiro-Wilk test (n_nodes=100, n_exam=200)
cc_ergraph: p=0.09412982314825058
cpl_ergraph: p=0.5377574563026428

Shapiro-Wilk test (n_nodes=100, n_exam=300)
cc_ergraph: p=0.418567955493927
cpl_ergraph: p=0.4255199432373047

Shapiro-Wilk test (n_nodes=100, n_exam=500)
cc_ergraph: p=0.4211519658565521
cpl_ergraph: p=0.8973879814147949
-----------------------------------------------

　それぞれでCrand、Lrandがいくらになったのかを求めておくのを忘れてしまったのだが、分布がどうであれCrand、Lrand自体はn_examが小さくてもある程度の狭い範囲に収まっているようにみえる。なので最初にざっくりデータの傾向をみる時点ではそこまで大きな回数のERグラフを生成する必要はないようにも思う。2通りしかやってないので雰囲気でしか話せないが、ノード数と同じくらいでいいかな。厳密にデータ解析しようとしたら、最低でもノード数の2倍～4,5倍ってとこですか。
　その辺も含めて結局のところ、どの程度のブレ幅までなら許容できるかをその都度考えなきゃダメってことですかね。大げさな話、解析したい元のネットワークたちのCnetが0.3だったり0.7だったりとぶれてたら、Crandが0.45だろうが0.55だろうが解析はできそうですし。逆に0.49と0.51の違いを扱いたい場合には、かなりの回数が必要そう。

　それと、ここに載せてない試行の中で、n_examが大きいにも関わらずシャピロウィルク検定で正規性がないとされたことがままあった。そういうときでもヒストグラム的にはとんがりコーンになっていたので、今回の目的に対してシャピロウィルク検定はいまいちよくないのかもしれない。てっきり回数稼げば正規分布になるだろうと思い込んでいたので、これについては近日中に考えようと思う。

　また、もしかしてCrand、Lrandはノード数の影響を受けないかも？とも思い始めた。や、ノード数2つしか試してないのでとやかく言えないのだが、雰囲気で話しているのでお許し願いたい。上ではリンク率0.5で通してたので、こちらもいろいろ変えてみる必要がありそう。

　今後の課題
・ノード数と必要なERグラフの数の関係性
・上のリンク率による影響
・リンク率を変えたときのCrand、Lrandの変化をみる
・weighted networkへの応用(もともとweighted networkを扱いたかったのだが、weightedなランダムネットワークをどのように作るかを考えた結果、まずunweightedな場合についてちゃんと知ってからweightedな解析について検証しよう、というのが今回のモチベである)