相関ネットワークを比較する

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

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

| 項番 | ページ内リンク |
|:----:|:-------------|
| 1 | [1. フィッシャー変換を用いた検定](#1-フィッシャー変換を用いた検定) |
| 2 | [2. パーミュテーション検定](#2-パーミュテーション検定) |


## 1. フィッシャー変換を用いた検定
フィッシャー変換は次式のように表すことができます。
```math
\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}
```
ここで、母相関係数ρについても同様に、
```math
\begin{flalign}
  &\eta = \frac{1}{2} \ln \frac{1 + \rho}{1 - \rho}&
\end{flalign}
```
とすることで、$\sqrt{n-3} (z-\eta)$ が正規分布に従うことが知られています。


また、各条件での相関係数の比較は、以下の式で求められるZスコアを用いて検定し、*p*値を算出することができます。
```math
\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で実装するとこんな感じでしょうか。
```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
```
実装には以下の記事を参考にさせて頂きました。

https://github.com/c60evaporator/covariance-visualization

https://toukeier.hatenablog.com/entry/correlateion-coefficient-difference-test

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

https://qiita.com/WolfMoon/items/ac217abd97483a5e3d49

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

以下が実装になります。
```python
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](https://github.com/groovy-phazuma/ML_DL_Notebook/tree/main/Correlation_Comparison)で公開しています。

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

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

https://www.coronasha.co.jp/np/isbn/9784339027327/