偏相関係数
疑似相関
2つの変数XとYの相関係数が、第３の変数Zの影響を除くと相関が無くなることがある。この疑似相関（spurious correlation）に対して、第３の変数Zの影響を単回帰モデルによって変数XとYから除き、その残りの部分について相関係数を求めたものは偏相関係数（partial correlation coefficient）と呼ばれている。すなわち、２つの単回帰モデル
において、
と
の相関係数を偏相関係数と呼び、
で表す。
である。
XとY、YとZ、XとZの間の相関係数を
、
、
で表すと次式が成り立つ。
＜岡本安晴「Delphiで学ぶデータ分析法」CQ出版社＞のpp．32-35に証明がある。
式（４）を導く考え方については、＜岡本安晴「データを使いこなすための統計入門」22世紀アート、電子書籍＞のp.40-41でも説明している。
ここでは、式（４）により代数的に偏相関係数を求めるのではなく、ベイズ分析により偏相関係数の事後分布を求める。事後分布から、疑似相関係数がどの範囲にあるかの推測を行うことができる。
式（１）、（２）、（３）に基づく偏相関係数を求めるStanスクリプトをリスト１のように用意した。リスト１では、変数XおよびYから変数Zの影響を除くのに単回帰モデル（式（１）および（２））を用いているが、これは一般化回帰モデルを設定することもできる。また、影響を除く変数も変数Zの１つではなく、複数の変数を設定することも可能である。
リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプトをリスト２のように用意した。
リスト２のプログラムでは、3変量の相関係数をリスト３のStanスクリプトによって求めている。リスト１では相関係数の事前分布として区間［−1，1］上の一様分布
を設定しているが、リスト３では相関係数の確率分布としてLKJ相関行列分布を用いている。リスト１では、２変数の間の相関係数なので一様分布として
が使えるが、リスト３では３変数の間の相関係数が扱われている。相関行列の成分としての相関係数には、相関行列の行列式の値が非負という制約がある。この制約のため、各相関係数に一様分布を独立に設定することはできない。３変数の場合に、３つの相関係数の値を座標とする３次元空間において、相関行列の行列式が正の領域と非正の領域を示す３D散布図を描くプログラムをこのウェブサイトに示した。LKJ相関分布は、パラメータ
のとき一様分布であるが、周辺分布は3次以上の相関行列では一様分布ではない。３次の相関行列の場合の周辺分布をこのウェブサイトで示した。同時分布が一様分布であっても周辺分布が一様分布にならない２次元分布の例もこのウェブサイトの後半で示している。３次の相関行列のLKJ分布の３次元散布図を描くプログラムはこのウェブサイトに用意した。
リスト１、リスト２，リスト３のファイルとサンプルデータファイル（図１、図２のもの）を圧縮ファイルSpurCorrFiles.zipにまとめた。ダウンロードして自由に利用されたい。
偏相関係数を求めるプログラムの実行例を以下に示すが、データは図１と図２の2通りを用意した。
図１ DataXYZ.csv
図２ DataXYZW.csv
図１、図２のデータはともにCSV形式のファイルである。CSV形式のファイルは、Excelで図１、図２のようにデータを用意して、保存時にファイル拡張子として「.csv」を選ぶとCSV形式のファイルとして保存される。データは１行目に変数名を書き、2行目からデータ値を書く。１列目はデータの識別用で計算には用いられないが、入力データの形式として必要である。図１のデータは、第３の変数Zの影響を変数XとYから除くと相関係数（偏相関係数）は０であると考えられる値になる。図２のデータは、変数XとYの相関係数をそのまま算出すると正の値であるが、第３の変数Zの効果を除くと、隠されていた負の相関関係が現れ、偏相関係数は負の値になる。
まず、図１のデータDataXYZ.csvを分析する。
リスト２のプログラムcalc_PartialCorr.pyを以下のように実行する。
(py39) PS D:\SpurCorrFiles> python
calcPCorr.py
INFO:numexpr.utils:NumExpr defaulting
to 8 threads.
Input file name (*.csv) = DataXYZ.csv
入力データファイル名を聞いてくるので、図１のデータファイルDataXYZ.csvを入力する。Enterキーを押すと計算が始まる。
計算が終了すると図３のグラフが表示される。
図３
グラフはファイル名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が出力されている。
XとYの相関係数は正（MAP推定値0.488；95% CI = [0.371, 0.585]）であるが、Zの影響を単回帰モデルで調整した偏相関係数は０（MAP推定値-0.032；95% CI=[-0.172, 0.104]）であると見なせる値である。
図２のデータDataXYZW.csvを分析すると、図４の結果を得る。
図４
出力ファイルに出力されている相関係数などは以下のようである。
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]
正であるXとYの相関係数（MAP推定値0.628；95% CI=[0.525, 0.699]）の偏相関係数は負（MAP推定値-0.518；95% CI=[-0.609, -0.404]）になっている。
リスト１ 偏相関係数の計算用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]]);
}
}
リスト２ 偏相関係数を求める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()
リスト３ ３変量の相関係数を求める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のインストールはこのウェブサイトで説明している。
Stan（PyStan）の説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。