はじめに

ネットワークはノードとエッジの集合として表現されます。特に、ノード間の相関係数をエッジとして用いる「相関ネットワーク」は幅広い分野で使用されます。
現実的なタスクとして、条件Aと条件Bの相関ネットワークを構築し比較することが考えられます。この際に、どのように相関ネットワークを比較すれば良いのでしょうか？

今回は、1. フィッシャー変換を用いた手法と2. パーミュテーション検定を用いた相関係数の差の明示的な検定について紹介します。

項番	ページ内リンク
1	1. フィッシャー変換を用いた検定
2	2. パーミュテーション検定

1. フィッシャー変換を用いた検定

フィッシャー変換は次式のように表すことができます。

\begin{flalign}
  &Z_A(x_i,x_j) = \frac{1}{2} \ln \frac{1 + Cor_A (x_i,x_j)}{1 - Cor_A (x_i,x_j)}&
\end{flalign}

ここで、母相関係数ρについても同様に、

\begin{flalign}
  &\eta = \frac{1}{2} \ln \frac{1 + \rho}{1 - \rho}&
\end{flalign}

とすることで、$\sqrt{n-3} (z-\eta)$ が正規分布に従うことが知られています。

また、各条件での相関係数の比較は、以下の式で求められるZスコアを用いて検定し、p値を算出することができます。

\begin{flalign}
  &Z_{AB}(x_i,x_j) = \frac{Z_A(x_i,x_j) - Z_B(x_i,x_j)}{\sqrt{\frac{1}{n_A - 3} + \frac{1}{n_B - 3}}}&
\end{flalign}

Pythonで実装するとこんな感じでしょうか。

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.datasets import load_iris

def Fisher_transformation(df=None,x='sepal length (cm)',y='petal length (cm)',p=0.05,verbose=True):
    if df is None:
        data = load_iris()
        df = pd.DataFrame(data['data'],columns=data['feature_names'])
    
    n = len(df)
    r = np.corrcoef(df[x], df[y])[0,1]
    z = np.log((1 + r) / (1 - r)) / 2

    eta_min = z - stats.norm.ppf(q=1-p/2, loc=0, scale=1) / np.sqrt(n - 3)
    eta_max = z - stats.norm.ppf(q=p/2, loc=0, scale=1) / np.sqrt(n - 3)

    rho_min = (np.exp(2 * eta_min) - 1) / (np.exp(2 * eta_min) + 1)
    rho_max = (np.exp(2 * eta_max) - 1) / (np.exp(2 * eta_max) + 1)

    if verbose:
        print(r)
        print(f'95% confident interval: {rho_min}〜{rho_max}')

    return z

def Cor_diff_test(df1,df2,x='sepal length (cm)',y='petal length (cm)',verbose=True):
    z1 = Fisher_transformation(df1,x=x,y=y,verbose=verbose)
    z2 = Fisher_transformation(df2,x=x,y=y,verbose=verbose)

    T = (z1 - z2) / np.sqrt((1/(len(df1)-3)) + (1/(len(df2) - 3)))
    p_value = stats.norm.cdf(T,loc=0,scale=1)

    return T, p_value

実装には以下の記事を参考にさせて頂きました。

2. パーミュテーション検定

さて、上記のp値の推定は、$Z_{AB}(x_i,x_j)$ が正規分布に従うことを仮定していますが、場合によっては成立しません。この場合はpermutation testを行い、経験的p値を算出する手法が有用です。
以下の記事が参考になるかもしれません。

ランダムに選択したサンプルをデータセット間で入れ替える操作を$N_{perm}$回実行し、各試行について取得した統計量と、元の順序での統計量の順位に着目することで、分布に依存せずにp値を算出することができます。

以下が実装になります。

import numpy as np
import pandas as pd
from tqdm import tqdm

def permutation_test(df1,df2,x='sepal length (cm)',y='petal length (cm)',n_perm=1000,alternative="less"):
    if alternative in ['less','greater']:
        pass
    else:
        raise ValueError("!! Inappropriate alternative type !!")
    original_t,original_p = ft.Cor_diff_test(df1,df2,x=x,y=y,verbose=False)

    n1 = len(df1)
    concat_df = pd.concat([df1,df2])

    # permutation
    perm_res = [original_t]
    for i in tqdm(range(n_perm)):
        shuffle_df = concat_df.sample(frac=1, random_state=i)
        u_df1 = shuffle_df.iloc[0:n1,:]
        u_df2 = shuffle_df.iloc[n1:,:]
        ut,up = ft.Cor_diff_test(u_df1,u_df2,x=x,y=y,verbose=False)
        perm_res.append(ut)
    
    # calc p-value
    if alternative == "less":
        perm_res = sorted(perm_res)
    else:
         perm_res = sorted(perm_res,reverse=True)

    original_idx = perm_res.index(original_t)
    perm_p = original_idx / n_perm

    return perm_p

おわりに

今回は相関ネットワークを比較する検定手法を紹介しました。コードはGitHubで公開しています。

本記事は竹本先生が書かれた「生物ネットワーク解析」を拝読し、大いに参考にさせていただきました [1]。相関ネットワーク解析に限らず、幅広い内容が収録されているため、興味がある方はぜひ一度お読みいただくことをおすすめします。

References

[1] 竹本和広 (2021). 生物ネットワーク解析　コロナ社