Up

偏相関係数

疑似相関

 

PyStan 3 対応版をこのウェブサイトに用意した。

 

 

2つの変数XYの相関係数が、第3の変数Zの影響を除くと相関が無くなることがある。この疑似相関(spurious correlation)に対して、第3の変数Zの影響を単回帰モデルによって変数XYから除き、その残りの部分について相関係数を求めたものは偏相関係数(partial correlation coefficient)と呼ばれている。すなわち、2つの単回帰モデル

において、の相関係数を偏相関係数と呼び、で表す。

である。

XYYZXZの間の相関係数をで表すと次式が成り立つ。

岡本安晴「Delphiで学ぶデータ分析法」CQ出版社>のpp32-35に証明がある。

式(4)を導く考え方については、<岡本安晴「データを使いこなすための統計入門」22世紀アート、電子書籍>のp.40-41でも説明している。

 

ここでは、式(4)により代数的に偏相関係数を求めるのではなく、ベイズ分析により偏相関係数の事後分布を求める。事後分布から、疑似相関係数がどの範囲にあるかの推測を行うことができる。

式(1)、(2)、(3)に基づく偏相関係数を求めるStanスクリプトをリスト1のように用意した。リスト1では、変数XおよびYから変数Zの影響を除くのに単回帰モデル(式(1)および(2))を用いているが、これは一般化回帰モデルを設定することもできる。また、影響を除く変数も変数Zの1つではなく、複数の変数を設定することも可能である。

リスト1のStanスクリプトを用いて偏相関係数を求めるPythonスクリプトをリスト2のように用意した。

リスト2のプログラムでは、3変量の相関係数をリスト3Stanスクリプトによって求めている。リスト1では相関係数の事前分布として区間[−11]上の一様分布を設定しているが、リスト3では相関係数の確率分布としてLKJ相関行列分布を用いている。リスト1では、2変数の間の相関係数なので一様分布としてが使えるが、リスト3では3変数の間の相関係数が扱われている。相関行列の成分としての相関係数には、相関行列の行列式の値が非負という制約がある。この制約のため、各相関係数に一様分布を独立に設定することはできない。3変数の場合に、3つの相関係数の値を座標とする3次元空間において、相関行列の行列式が正の領域と非正の領域を示す3D散布図を描くプログラムをこのウェブサイトに示した。LKJ相関分布は、パラメータのとき一様分布であるが、周辺分布は3次以上の相関行列では一様分布ではない。3次の相関行列の場合の周辺分布をこのウェブサイトで示した。同時分布が一様分布であっても周辺分布が一様分布にならない2次元分布の例もこのウェブサイトの後半で示している。3次の相関行列のLKJ分布の3次元散布図を描くプログラムはこのウェブサイトに用意した。

リスト1、リスト2,リスト3のファイルとサンプルデータファイル(図1、図2のもの)を圧縮ファイルSpurCorrFiles.zipにまとめた。ダウンロードして自由に利用されたい。

 

偏相関係数を求めるプログラムの実行例を以下に示すが、データは図1と図2の2通りを用意した。

図1 DataXYZ.csv

 

図2 DataXYZW.csv

 

図1、図2のデータはともにCSV形式のファイルである。CSV形式のファイルは、Excelで図1、図2のようにデータを用意して、保存時にファイル拡張子として「.csv」を選ぶとCSV形式のファイルとして保存される。データは1行目に変数名を書き、2行目からデータ値を書く。1列目はデータの識別用で計算には用いられないが、入力データの形式として必要である。図1のデータは、第3の変数Zの影響を変数XYから除くと相関係数(偏相関係数)は0であると考えられる値になる。図2のデータは、変数XYの相関係数をそのまま算出すると正の値であるが、第3の変数Zの効果を除くと、隠されていた負の相関関係が現れ、偏相関係数は負の値になる。

まず、図1のデータDataXYZ.csvを分析する。

リスト2のプログラムcalc_PartialCorr.pyを以下のように実行する。

 

(py39) PS D:\SpurCorrFiles> python calcPCorr.py

INFO:numexpr.utils:NumExpr defaulting to 8 threads.

Input file name (*.csv) = DataXYZ.csv

 

入力データファイル名を聞いてくるので、図1のデータファイルDataXYZ.csvを入力する。Enterキーを押すと計算が始まる。

計算が終了すると図3のグラフが表示される。

図3

 

グラフはファイル名Figure.pngで保存される。計算結果の出力ファイルはResults.txtである。データDataXYZ.csvを入力ファイルとした場合、出力ファイルの最後に相関係数などが以下のように出力されている。

 

Correlation Coefficient Rxy

  MAP Est. = 0.488

  Mean     = 0.485

  Median   = 0.487

  95% CI   = [0.371, 0.585]

 

Partial Correlation Coefficient Rxy.z

  MAP Est. = -0.032

  Mean     = -0.036

  Median   = -0.035

  95% CI   = [-0.172, 0.104]

 

Correlation Coefficient Rxz

  MAP Est. = 0.712

  Mean     = 0.705

  Median   = 0.707

  95% CI   = [0.629, 0.769]

 

Correlation Coefficient Ryz

  MAP Est. = 0.720

  Mean     = 0.712

  Median   = 0.714

  95% CI   = [0.637, 0.775]

 

MAP推定値、平均値、中央値および95 CIが出力されている。

XYの相関係数は正(MAP推定値0.48895% CI = [0.371, 0.585])であるが、Zの影響を単回帰モデルで調整した偏相関係数は0(MAP推定値-0.03295% CI=[-0.172, 0.104])であると見なせる値である。

 

図2のデータDataXYZW.csvを分析すると、図4の結果を得る。

図4

 

出力ファイルに出力されている相関係数などは以下のようである。

 

Correlation Coefficient Rxy

  MAP Est. = 0.628

  Mean     = 0.618

  Median   = 0.621

  95% CI   = [0.525, 0.699]

 

Partial Correlation Coefficient Rxy.z

  MAP Est. = -0.518

  Mean     = -0.511

  Median   = -0.513

  95% CI   = [-0.609, -0.404]

 

Correlation Coefficient Rxz

  MAP Est. = 0.868

  Mean     = 0.865

  Median   = 0.866

  95% CI   = [0.827, 0.897]

 

Correlation Coefficient Ryz

  MAP Est. = 0.868

  Mean     = 0.864

  Median   = 0.865

  95% CI   = [0.825, 0.897]

 

正であるXYの相関係数(MAP推定値0.62895% CI=[0.525, 0.699])の偏相関係数は負(MAP推定値-0.51895% CI=[-0.609, -0.404])になっている。

 

 

リスト1 偏相関係数の計算用Stanスクリプト(ファイル名PartialCorr.stan

 

data {

    int N;

    real X[N];

    real Y[N];

    real Z[N];

}

parameters {

    real mu1;

    real mu2;

    real<lower = 0.0> sgm1;

    real<lower = 0.0> sgm2;

    real<lower = 0.0> sgm_x;

    real<lower = 0.0> sgm_y;

    real<lower = -1.0, upper = 1.0> rho;

    real Ax;

    real Bx;

    real Ay;

    real By;

}

transformed parameters {

    vector[2] XY[N];

    for (i in 1:N) {

        XY[i][1] = X[i] - (Bx + Ax * Z[i]);

        XY[i][2] = Y[i] - (By + Ay * Z[i]);  

    }

}

model {

    mu1 ~ normal(0.0, 100.0);

    mu2 ~ normal(0.0, 100.0);

    sgm1 ~ exponential(0.001);

    sgm2 ~ exponential(0.001);

    rho ~ uniform(-1.0, 1.0);

    Ax ~ normal(0.0, 1000.0);

    Ay ~ normal(0.0, 1000.0);

    Bx ~ normal(0.0, 1000.0);

    By ~ normal(0.0, 1000.0);

    sgm_x ~ exponential(0.001);

    sgm_y ~ exponential(0.001);

    for (i in 1:N) {

        X[i] ~ normal(Bx + Ax * Z[i], sgm_x);

        Y[i] ~ normal(By + Ay * Z[i], sgm_y);

        XY[i] ~ multi_normal([mu1, mu2], [[sgm1*sgm1, rho*sgm1*sgm2], [rho*sgm2*sgm1, sgm2*sgm2]]);

    }

}

 

 

リスト2 偏相関係数を求めるPythonスクリプト(ファイル名calcPCorr.py

 

import csv

import numpy as np

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

import scipy.stats as ss

import pickle

 

fin_nm = input('Input file name (*.csv) = ')

with open(fin_nm, 'r') as f:

    csv_data = [d for d in csv.reader(f)]

 

savefig_nm = 'Figure.png'

fout_nm = 'Results.txt'

fout = open(fout_nm, 'w')

fout.write('Input data file = {}'.format(fin_nm))

 

N = len(csv_data) - 1

print('N =', N)

X = []

Y = []

Z = []

for i in range(N):

    X.append(float(csv_data[i+1][1]))

    Y.append(float(csv_data[i+1][2]))

    Z.append(float(csv_data[i+1][3]))

 

fout.write('\n\nData =\n')

for x, y, z in zip(X, Y, Z):

    fout.write(f'{x}  {y}  {z}\n')

 

print('Compiling the stan script for the correlation coefficients.')

sm0 = pystan.StanModel(file = 'LKJCorr.stan')

 

Data = {'N': len(X), 'X':X, 'Y':Y, 'Z':Z}

 

print('MCMC sampling for the correlations started.')

 

 

import platform

if platform.system() == 'Windows':

    fit0 = sm0.sampling(data = Data,

                        pars = ['muX', 'muY', 'muZ', 'sgmX', 'sgmY', 'sgmZ', 'RhoXY', 'RhoXZ', 'RhoYZ'],

                        iter = 5000, n_jobs = 1)

else:

    fit0 = sm0.sampling(data = Data,

                        pars = ['muX', 'muY', 'muZ', 'sgmX', 'sgmY', 'sgmZ', 'RhoXY', 'RhoXZ', 'RhoYZ'],

                        iter = 5000)

 

print(fit0)

fout.write('\nMCMC sampling summary (Three variate normal distribution)\n')

fout.write('{}\n'.format(fit0))

 

 

print('Compiling the stan script for the partial correlation coefficient.')

sm1 = pystan.StanModel(file = 'PartialCorr.stan')

 

print('MCMC sampling for the partial correlation coefficient started.')

if platform.system() == 'Windows':

    fit1 = sm1.sampling(data = Data, pars = ['mu1', 'mu2', 'sgm1', 'sgm2', 'sgm_x', 'sgm_y',

                                      'rho', 'Ax', 'Bx', 'Ay', 'By'], iter = 5000, n_jobs = 1)    

else:

    fit1 = sm1.sampling(data = Data, pars = ['mu1', 'mu2', 'sgm1', 'sgm2', 'sgm_x', 'sgm_y',

                                      'rho', 'Ax', 'Bx', 'Ay', 'By'], iter = 5000)    

 

 

print(fit1)

fout.write('\nMCMC sampling summary (Partial correlation coefficient)\n')

fout.write('{}\n'.format(fit1))

 

 

RhoXY = fit0['RhoXY']

RhoXZ = fit0['RhoXZ']

RhoYZ = fit0['RhoYZ']

p_corr = fit1['rho']

 

 

def CalcMAPEst(samples, a = 0.05, n_points = 10000):

    """

        Calculatte a MAP estimate from a KDE graph on [Lp, Up]

        Lp and Up are 100*a/2 and 100(1-a/2) percentile points of samples

    """

    Lp, Up = np.percentile(samples, [100 * a/2, 100 * (1 - a/2)])  # import numpy as np

    coord = np.linspace(Lp, Up, n_points) 

    est_pdf = ss.gaussian_kde(samples).pdf(coord)  # import scipy.stats as ss

    map_idx = np.argmax(est_pdf) 

    MAP_Est = coord[map_idx]       

    return MAP_Est, est_pdf[map_idx]

 

#       日本語表記の設定

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

 

plt.subplots(figsize = (9, 6))

 

xcoord = np.linspace(-1, 1, 1000)

ycoord = ss.gaussian_kde(RhoXY).pdf(xcoord)

r_xy_map = CalcMAPEst(RhoXY)[0]

r_xy_mean = RhoXY.mean()

r_xy_025, r_xy_med, r_xy_975 = np.percentile(RhoXY, [2.5, 50, 97.5])

print('\n相関係数 Rxy')

print('MAP推定値 = {0:.3f}'.format(r_xy_map))

print('   平均値 = {0:.3f}'.format(r_xy_mean))

print('   中央値 = {0:.3f}'.format(r_xy_med))

print('   95% CI = [{0:.3f}, {1:.3f}]'.format(r_xy_025, r_xy_975))

 

fout.write('\nCorrelation Coefficient Rxy\n')

fout.write('  MAP Est. = {0:.3f}\n'.format(r_xy_map))

fout.write('  Mean     = {0:.3f}\n'.format(r_xy_mean))

fout.write('  Median   = {0:.3f}\n'.format(r_xy_med))

fout.write('  95% CI   = [{0:.3f}, {1:.3f}]\n'.format(r_xy_025, r_xy_975))

 

plt.subplot(2,2,1)

plt.title(('相関係数 ' + r'$\rho_{xy}$') + (': mode={0:.3f}'.format(r_xy_map)),

          fontsize = 16)

plt.plot(xcoord, ycoord)

 

ycoord = ss.gaussian_kde(p_corr).pdf(xcoord)

p_corr_map = CalcMAPEst(p_corr)[0]

p_corr_mean = p_corr.mean()

p_corr_025, p_corr_med, p_corr_975 = np.percentile(p_corr, [2.5, 50, 97.5])

print('\n偏相関係数 Rxy.z')

print('MAP推定値 = {0:.3f}'.format(p_corr_map))

print('   平均値 = {0:.3f}'.format(p_corr_mean))

print('   中央値 = {0:.3f}'.format(p_corr_med))

print('   95% CI = [{0:.3f}, {1:.3f}]'.format(p_corr_025, p_corr_975))

 

fout.write('\nPartial Correlation Coefficient Rxy.z\n')

fout.write('  MAP Est. = {0:.3f}\n'.format(p_corr_map))

fout.write('  Mean     = {0:.3f}\n'.format(p_corr_mean))

fout.write('  Median   = {0:.3f}\n'.format(p_corr_med))

fout.write('  95% CI   = [{0:.3f}, {1:.3f}]\n'.format(p_corr_025, p_corr_975))

 

plt.subplot(2,2,2)

plt.title('偏相関係数 ' + r'$\rho_{XY{\bullet}Z}$' +

          ': mode={0:.3f}'.format(p_corr_map), fontsize = 16)

plt.plot(xcoord, ycoord)

 

ycoord = ss.gaussian_kde(RhoXZ).pdf(xcoord)

r_xz_map = CalcMAPEst(RhoXZ)[0]

r_xz_mean = RhoXZ.mean()

r_xz_025, r_xz_med, r_xz_975 = np.percentile(RhoXZ, [2.5, 50, 97.5])

print('\n相関係数 Rxz')

print('MAP推定値 = {0:.3f}'.format(r_xz_map))

print('   平均値 = {0:.3f}'.format(r_xz_mean))

print('   中央値 = {0:.3f}'.format(r_xz_med))

print('   95% CI = [{0:.3f}, {1:.3f}]'.format(r_xz_025, r_xz_975))

 

fout.write('\nCorrelation Coefficient Rxz\n')

fout.write('  MAP Est. = {0:.3f}\n'.format(r_xz_map))

fout.write('  Mean     = {0:.3f}\n'.format(r_xz_mean))

fout.write('  Median   = {0:.3f}\n'.format(r_xz_med))

fout.write('  95% CI   = [{0:.3f}, {1:.3f}]\n'.format(r_xz_025, r_xz_975))

 

plt.subplot(2,2,3)

plt.title(r'$\rho_{xz}$' + ': mode={0:.3f}'.format(r_xz_map),

          fontsize = 16)

plt.plot(xcoord, ycoord)

 

ycoord = ss.gaussian_kde(RhoYZ).pdf(xcoord)

r_yz_map = CalcMAPEst(RhoYZ)[0]

r_yz_mean = RhoYZ.mean()

r_yz_025, r_yz_med, r_yz_975 = np.percentile(RhoYZ, [2.5, 50, 97.5])

print('\n相関係数 Ryz')

print('MAP推定値 = {0:.3f}'.format(r_yz_map))

print('   平均値 = {0:.3f}'.format(r_yz_mean))

print('   中央値 = {0:.3f}'.format(r_yz_med))

print('   95% CI = [{0:.3f}, {1:.3f}]'.format(r_yz_025, r_yz_975))

 

fout.write('\nCorrelation Coefficient Ryz\n')

fout.write('  MAP Est. = {0:.3f}\n'.format(r_yz_map))

fout.write('  Mean     = {0:.3f}\n'.format(r_yz_mean))

fout.write('  Median   = {0:.3f}\n'.format(r_yz_med))

fout.write('  95% CI   = [{0:.3f}, {1:.3f}]\n'.format(r_yz_025, r_yz_975))

 

plt.subplot(2,2,4)

plt.title(r'$\rho_{yz}$' + ': mode={0:.3f}'.format(r_yz_map),

          fontsize = 16)

plt.plot(xcoord, ycoord)

 

plt.tight_layout()

 

plt.savefig(savefig_nm)

print('\nFile {} of the figure was saved.'.format(savefig_nm))

plt.show()

 

fout.close()

print('\nFile {} was saved.'.format(fout_nm))

print()

 

 

リスト3 3変量の相関係数を求めるStanスクリプト(ファイル名LKJCorr.stan

 

data {

    int N;

    real X[N];

    real Y[N];

    real Z[N];

}

parameters {

    real muX;

    real muY;

    real muZ;

    real<lower = 0.0> sgmX;

    real<lower = 0.0> sgmY;

    real<lower = 0.0> sgmZ;

    corr_matrix[3] cor;

}

transformed parameters {

    real RhoXY;

    real RhoXZ;

    real RhoYZ;

    cov_matrix[3] Cov;

    Cov[1][1] = sgmX*sgmX;

    Cov[1][2] = cor[1][2]*sgmX*sgmY;

    Cov[1][3] = cor[1][3]*sgmX*sgmZ;

    Cov[2][1] = cor[2][1]*sgmY*sgmX;

    Cov[2][2] = sgmY*sgmY;

    Cov[2][3] = cor[2][3]*sgmY*sgmZ;

    Cov[3][1] = cor[3][1]*sgmX*sgmZ;

    Cov[3][2] = cor[3][2]*sgmY*sgmZ;

    Cov[3][3] = sgmZ*sgmZ;

    RhoXY = cor[1][2];

    RhoXZ = cor[1][3];

    RhoYZ = cor[2][3];

}

model {

    muX ~ normal(0.0, 1000.0);

    muY ~ normal(0.0, 1000.0);

    muZ ~ normal(0.0, 1000.0);

    sgmX ~ exponential(0.001);

    sgmY ~ exponential(0.001);

    sgmZ ~ exponential(0.001);

    cor ~ lkj_corr(1.0);

    for (i in 1:N) {

        [X[i], Y[i], Z[i]] ~ multi_normal([muX, muY, muZ], Cov);

    }

}

 

 

PyStanのインストールはこのウェブサイトで説明している。

StanPyStan)の説明は、<岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版>で行っている。

 

 

Up