偏相関係数（PyStan 3対応）
疑似相関
サンプルスクリプトは、
Ubuntu 22/WSL/Windows 10において、Python 3.11/Anaconda環境で
Pystan 3.10.0、Arviz 0.18.0を用いて実行した。
2つの変数XとYの相関係数が、第３の変数Zの影響を除くと相関が無くなることがある。この疑似相関（spurious correlation）に対して、第３の変数Zの影響を単回帰モデルによって変数XとYから除き、その残りの部分について相関係数を求めたものは偏相関係数（partial correlation coefficient）と呼ばれている。すなわち、２つの単回帰モデル
において、
と
の相関係数を偏相関係数と呼び、
で表す。
である。
XとY、YとZ、XとZの間の相関係数を
、
、
で表すと次式が成り立つ。
＜岡本安晴「Delphiで学ぶデータ分析法」CQ出版社＞のpp．32-35に証明がある。
ここでは、式（４）により代数的に偏相関係数を求めるのではなく、ベイズ分析により偏相関係数の事後分布を求める。事後分布から、偏相関係数の存在範囲の推測を行うことができる。
式（１）、（２）、（３）に基づく偏相関係数を求めるStanスクリプトをリスト１のように用意した。リスト１では、変数XおよびYから変数Zの影響を除くのに単回帰モデル（式（１）および（２））を用いている。すなわち、
for (i in 1:N) {
XY[i][1] = X[i] - (Bx + Ax * Z[i]); // Residual
XY[i][2]
= Y[i] - (By + Ay * Z[i]);
// Residual
}
・・・
for (i in 1:N) {
X[i]
~ normal(Bx + Ax * Z[i], sgm_x);
Y[i]
~ normal(By + Ay * Z[i], sgm_y);
//
// The
model for estimation of the partial correlation rho
//
XY[i] ~ multi_normal([mu1, mu2], [[sgm1*sgm1, rho*sgm1*sgm2],
[rho*sgm2*sgm1, sgm2*sgm2]]);
である。
これは一般化回帰モデルを設定することもできる。また、影響を除く変数も変数Zの１つではなく、複数の変数を設定することも可能である。
リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプトをリスト２のように用意した。
リスト２のプログラムでは、3変量の相関係数をリスト３のStanスクリプトによって求めている。
リスト１、リスト２，リスト３のファイルとサンプルデータファイル（図１、図２のもの）を圧縮ファイルSpurCorrPS3Files.zipにまとめた。ダウンロードして自由に利用されたい。
偏相関係数を求めるプログラムの実行例を以下に示す。サンプルスクリプトは、Ubuntu 22/WSL/Windows 10において、Python 3.11/Anaconda環境でPystan 3.10.0、Arviz 0.18.0を用いて実行した。
データは図１と図２の2通りを用意した。
図１ DataXYZ.csv
図２ DataXYZW.csv
図１、図２のデータはともにCSV形式のファイルである。CSV形式のファイルは、Excelで図１、図２のようにデータを用意して、保存時にファイル拡張子として「.csv」を選ぶとCSV形式のファイルとして保存される。データは１行目に変数名を書き、2行目からデータ値を書く。１列目はデータの識別用で計算には用いられないが、入力データの形式として必要である。
図１のデータは、第３の変数Zの影響を変数XとYから除くと相関係数（偏相関係数）は０であると考えられる値になる。
図２のデータは、変数XとYの相関係数をそのまま算出すると正の値であるが、第３の変数Zの効果を除くと、隠されていた負の相関関係が現れ、偏相関係数は負の値になる。
まず、図１のデータDataXYZ.csvを分析する。
リスト２のプログラムcalcPCorr.pyを以下のように実行する。
(py311)
yasuharu@PSY-PM-PC:/*****/SpurCorrPS3Files$ python calcPCorr.py
Input file name (*.csv) = DataXYZ.csv
入力データファイル名を聞いてくるので、図１のデータファイルDataXYZ.csvを入力する。Enterキーを押すと計算が始まる。
計算が終了すると図３のグラフが表示される。
事後分布の点推定値として中央値が示されている。Gelmanら（2021）は、事後分布の点推定値として中央値を薦めている。点推定値としては、MAP推定値がよいという報告（岡本、2022）があるが、計算コストを考えて、ここでは中央値を用いた。
図３
グラフは、ファイル名Figure.pngで保存される。計算結果の出力ファイルは、Results.txtである。データDataXYZ.csvを入力ファイルとした場合、出力ファイルの最後に相関係数などが以下のように出力されている。
Correlation Coefficient Rxy
Mean =
0.488
Median = 0.489
95% CI = [0.373,
0.587]
Partial Correlation Coefficient Rxy.z
Mean =
-0.035
Median = -0.034
95% CI = [-0.168,
0.100]
Correlation Coefficient Rxz
Mean =
0.707
Median = 0.709
95% CI = [0.630,
0.771]
Correlation Coefficient Ryz
Mean =
0.714
Median = 0.715
95% CI = [0.640, 0.777]
平均値、中央値および95％ CIが出力されている。
XとYの相関係数は正（中央値0.489；95% CI = [0.373, 0.587]）であるが、Zの影響を単回帰モデルで調整した偏相関係数は０（中央値-0.035；95% CI=[-0.168, 0.100]）であると見なせる値である。
次に、図２のデータDataXYZW.csvを分析する。図４の結果を得る。
図４
出力ファイルに出力されている相関係数などは以下のようである。
Correlation Coefficient Rxy
Mean =
0.619
Median = 0.619
95% CI = [0.529,
0.700]
Partial Correlation Coefficient Rxy.z
Mean =
-0.510
Median = -0.511
95% CI = [-0.608,
-0.402]
Correlation Coefficient Rxz
Mean =
0.865
Median = 0.866
95% CI = [0.828,
0.897]
Correlation Coefficient Ryz
Mean =
0.864
Median = 0.865
95% CI = [0.828,
0.898]
正であるXとYの相関係数（中央値0.619；95% CI=[0.529, 0.700]）に対して、Zの影響を除いた偏相関係数は負（中央値-0.511；95% CI=[-0.608, -0.402]）になっている。
変数間の関係を分析するとき、第３の変数の影響の可能性があれば、その第３の変数の影響を考慮して分析する必要がある。
Gelman, A., Hill, J., and Vehtari, A. (2021) Regression and other stories. Cambridge University Press.
岡本安晴（1998）「Delphiで学ぶデータ分析法」CQ出版社
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
岡本安晴（2022）ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――、基礎心理学研究、第41巻、1-7．
リスト１ 偏相関係数の計算用Stanスクリプト（ファイル名PartialCorr.stan）
data {
int N;
array[N] real X;
array[N] real Y;
array[N] real Z;
}
parameters {
real vmu1;
real vmu2;
real<lower = 0.0>
vsgm1;
real<lower = 0.0>
vsgm2;
real<lower = 0.0>
vsgm_x;
real<lower = 0.0>
vsgm_y;
real<lower = -1.0, upper
= 1.0> rho;
real Ax;
real Bx;
real Ay;
real By;
}
transformed 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;
array[N] vector[2] XY;
mu1 = vmu1 * 100;
mu2 = vmu2 * 100;
sgm1 = 0.001 + vsgm1 * 1000;
sgm2 = 0.001 + vsgm2 * 1000;
sgm_x = 0.001 + vsgm_x *
1000;
sgm_y = 0.001 + vsgm_y *
1000;
for (i in 1:N) {
XY[i][1] = X[i] - (Bx + Ax * Z[i]); // Residual
XY[i][2] = Y[i] - (By + Ay * Z[i]); // Residual
}
}
model {
vmu1 ~ normal(0.0, 1.0);
vmu2 ~ normal(0.0, 1.0);
vsgm1 ~ exponential(0.1);
vsgm2 ~ exponential(0.1);
rho ~ uniform(-1.0,
1.0); // Partial correlation
Ax ~ normal(0.0, 10.0);
Ay ~ normal(0.0, 10.0);
Bx ~ normal(0.0, 10.0);
By ~ normal(0.0, 10.0);
vsgm_x ~ exponential(0.1);
vsgm_y ~ exponential(0.1);
for (i in 1:N) {
X[i]
~ normal(Bx + Ax * Z[i], sgm_x);
Y[i]
~ normal(By + Ay * Z[i], sgm_y);
//
//
The model for estimation of the partial correlation rho
//
XY[i] ~ multi_normal([mu1, mu2], [[sgm1*sgm1, rho*sgm1*sgm2],
[rho*sgm2*sgm1, sgm2*sgm2]]);
}
}
リスト２ 偏相関係数を求めるPythonスクリプト（ファイル名calcPCorr.py）
import csv
import numpy as np
import stan
import arviz as az
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
import japanize_matplotlib
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')
#
#
相関係数の計算
#
Data = {'N': len(X), 'X':X, 'Y':Y,
'Z':Z}
def f_init():
return {'vmuX':0.5,
'vmuY':0.5, 'vmuZ':0.5,
'vsgmX':0.5,
'vsgmY':0.5, 'vsgmZ':0.5,
'cl':[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] }
with open('LKJCorr.stan') as f_stan:
sm0 =
stan.build(f_stan.read(), data=Data)
fit0 = sm0.sample(num_chains=4,
init=[f_init()]*4)
print('\nMCMC sampling for the
correlations started.')
i_data0 =
az.from_pystan(posterior=fit0, posterior_model=sm0)
print(az.summary(i_data0))
fout.write('\nMCMC sampling summary
(Three variate normal distribution)\n')
fout.write('{}\n'.format(az.summary(i_data0)))
#
# 偏相関係数の計算
#
print('\nCompiling the stan script for
the partial correlation coefficient.')
with open('PartialCorr.stan') as
f_stan:
sm1 =
stan.build(f_stan.read(), data=Data)
fit1 = sm1.sample()
i_data1 =
az.from_pystan(posterior=fit1, posterior_model=sm1)
print(az.summary(i_data1))
fout.write('\nMCMC sampling summary
(Partial correlation coefficient)\n')
fout.write('{}\n'.format(az.summary(i_data1)))
#
# （偏）相関係数の出力とグラフ描画
#
RhoXY = fit0['RhoXY']
RhoXZ = fit0['RhoXZ']
RhoYZ = fit0['RhoYZ']
p_corr = fit1['rho']
plt.subplots(figsize = (9, 6))
xcoord = np.linspace(-1, 1, 1000)
ycoord =
ss.gaussian_kde(RhoXY).pdf(xcoord)
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(' 平均値 = {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(' 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}$') + (':
med.={0:.3f}'.format(r_xy_med)),
fontsize = 16)
plt.plot(xcoord, ycoord)
ycoord =
ss.gaussian_kde(p_corr).pdf(xcoord)
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(' 平均値 = {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(' 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}$' +
': med.={0:.3f}'.format(p_corr_med), fontsize = 16)
plt.plot(xcoord, ycoord)
ycoord =
ss.gaussian_kde(RhoXZ).pdf(xcoord)
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(' 平均値 = {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(' 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}$' + ':
med.={0:.3f}'.format(r_xz_med),
fontsize = 16)
plt.plot(xcoord, ycoord)
ycoord =
ss.gaussian_kde(RhoYZ).pdf(xcoord)
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(' 平均値 = {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(' 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}$' + ':
med.={0:.3f}'.format(r_yz_med),
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()
リスト３ ３変量の相関係数を求めるStanスクリプト（ファイル名LKJCorr.stan）
data {
int N;
array[N] real X;
array[N] real Y;
array[N] real Z;
}
parameters {
real vmuX;
real vmuY;
real vmuZ;
real<lower = 0.0>
vsgmX;
real<lower = 0.0>
vsgmY;
real<lower = 0.0>
vsgmZ;
cholesky_factor_corr[3] cL;
}
transformed parameters {
real muX;
real muY;
real muZ;
real RhoXY;
real RhoXZ;
real RhoYZ;
real<lower = 0.0>
sgmX;
real<lower = 0.0>
sgmY;
real<lower = 0.0>
sgmZ;
corr_matrix[3] Corr;
cov_matrix[3] Cov;
muX = vmuX * 100.0;
muY = vmuY * 100.0;
muZ = vmuZ * 100.0;
sgmX = 0.01 + vsgmX * 100;
sgmY = 0.01 + vsgmY * 100;
sgmZ = 0.01 + vsgmZ * 100;
//
// Correlation matrix
//
Corr = cL * cL';
//
// Covariance
matrix
//
Cov[1][1] = sgmX*sgmX;
Cov[1][2] =
Corr[1][2]*sgmX*sgmY;
Cov[1][3] =
Corr[1][3]*sgmX*sgmZ;
Cov[2][1] =
Corr[2][1]*sgmY*sgmX;
Cov[2][2] = sgmY*sgmY;
Cov[2][3] =
Corr[2][3]*sgmY*sgmZ;
Cov[3][1] =
Corr[3][1]*sgmX*sgmZ;
Cov[3][2] =
Corr[3][2]*sgmY*sgmZ;
Cov[3][3] = sgmZ*sgmZ;
//
//
Correlations
//
RhoXY = Corr[1][2];
RhoXZ = Corr[1][3];
RhoYZ = Corr[2][3];
}
model {
vmuX ~ normal(0.0, 10.0);
vmuY ~ normal(0.0, 10.0);
vmuZ ~ normal(0.0, 10.0);
vsgmX ~ exponential(0.1);
vsgmY ~ exponential(0.1);
vsgmZ ~ exponential(0.1);
cL ~ lkj_corr_cholesky(1.0);
for (i in 1:N) {
[X[i], Y[i], Z[i]] ~ multi_normal([muX, muY, muZ], Cov);
}
}