偏相関係数
疑似相関
CmdStanPy対応：2025.11
2つの変数XとYの相関係数が、第３の変数Zの影響を除くと相関が無くなることがある。この疑似相関（spurious correlation）に対して、第３の変数Zの影響を単回帰モデルによって変数XとYから除き、その残りの部分について相関係数を求めたものは偏相関係数（partial correlation coefficient）と呼ばれている（Mardia ei al., 1979）。すなわち、２つの単回帰モデル
において、
と
の相関係数を偏相関係数と呼び、
で表す。
である。
XとY、YとZ、XとZの間の相関係数を
、
、
で表すと次式が成り立つ。
＜岡本安晴「Delphiで学ぶデータ分析法」CQ出版社＞のpp．32-35に証明がある。
式（４）を導く考え方については、＜岡本安晴「データを使いこなすための統計入門」22世紀アート、電子書籍＞のp.40-41でも説明している。
ここでは、式（４）により代数的に偏相関係数を求めるのではなく、ベイズ分析により偏相関係数の事後分布を求める。事後分布から、偏相関係数がどの範囲にあるかの推測を行うことができる。
定義より、偏相関係数を
と置くとき、式（１）と式（２）に対して
と置くことができる。
上式に基づき、リスト1に示すStanスクリプトを用意した。
リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプトをリスト２に用意した。
ファイルはSpurCorrFiles.zipにまとめた。
偏相関係数を求めるプログラムの実行例を以下に示すが、データは図１と図２の2通りを用意した。
図１ DataXYZ.csv
図２ DataXYZW.csv
図１、図２のデータはともにCSV形式のファイルである。CSV形式のファイルは、Excelで図１、図２のようにデータを用意して、保存時にファイル拡張子として「.csv」を選ぶとCSV形式のファイルとして保存される。データは１行目に変数名を書き、2行目からデータ値を書く。１列目はデータの識別用で計算には用いられないが、入力データの形式として必要である。図１のデータは、第３の変数Zの影響を変数XとYから除くと相関係数（偏相関係数）は０であると考えられる値になる。図２のデータは、変数XとYの相関係数をそのまま算出すると正の値であるが、第３の変数Zの効果を除くと、隠されていた負の相関関係が現れ、偏相関係数は負の値になる。
まず、図１のデータDataXYZ.csvを分析する。
リスト１のStanスクリプトファイル（pcor.stan）、リスト２のPythonスクリプトファイル（main_calc.py）、および入力データファイ（DataXYZ.csv）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/SpurCorrFiles$ python main_calc.py
上のコマンドを実行すると、入力データファイル名の設定が求められる。
(stan) *****/SpurCorrFiles$
python main_calc.py
Input file name (*.csv) = DataXYZ.csv
上の例では、図１のファイル名が設定されている。
ファイル名が設定されると、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まり、終了すると、トレース図が表示される（図３）。
図３
図３のWindowを閉じると、図４のグラフが表示される。
図４
回帰直線のパラメータの事後分布のグラフである。
図４のWindowを閉じると、図５のグラフが表示される。
図５
偏相関係数の事後分布のグラフである。偏相関係数が正である事後確率は、0.30925であると表示されている。事後分布は、ほぼ原点を中心にして分布していると言える。
図５のWindowを閉じると、図６のグラフが表示される。
図６
点（Z,X）の散布図が左上、点（X.Z、Y.Z）の散布図が右上に示されている。X.ZとY.Zは、式（１）における残差と式（２）における残差を表している。下段の散布図は、点（X,Y）の散布図が左側、点（Z,Y）の散布図が右側に示されている。散布図における赤の実線は、回帰直線である。散布図の上に表示されているｒの値は相関係数である。
変量XとYの相関係数はr=0.496であるが、偏相関係数はr=0.036であり、図５の事後分布から、偏相関係数は０であると言えそうである。
図６のWindowを閉じると、スクリプトの実行終了である。
次に、図２のデータを分析してみる。
リスト２のスクリプトを実行して、入力データファイルを以下のように設定する。
(stan) *****/SpurCorrFiles$
python main_calc.py
Input file name (*.csv) = DataXYZW.csv
上記ファイル設定で、図７のトレース図が表示される。
図７
図７のWindowを閉じると、図８のグラフが表示される。
図８
図８のWindowを閉じると、図９のグラフが表示される。
図９
偏相関係数の事後分布のグラフであるが、偏相関係数が正である事後確率は、0.00000と表示されている。事後分布の中央値は、―0.511である。偏相関係数は、負であると言えそうである。
図９のWindowを閉じると、図10のグラフが表示される。
図10
左下の（X,Y）の散布図では、相関係数は0.629であり、回帰直線は右上がりの直線である。変量Zの影響を除いた点（X.Z、Y.Z）の散布図では、偏相関係数はー0.511で、回帰直線は右下がりの直線である。
図10のWindowを閉じると、スクリプトの実行終了である。
リスト１ 偏相関係数のStanスクリプト（pcor.stan）
data {
int N;
array[N] real X;
array[N] real Y;
array[N] real Z;
}
parameters {
real<lower = 0.001,
upper=1000> sgm1;
real<lower = 0.001,
upper=1000> sgm2;
real<lower = -0.999,
upper = 0.999> rho;
real Ax;
real Bx;
real Ay;
real By;
}
transformed parameters {
array[N] vector[2] XY;
for (i
in 1:N) {
XY[i][1] = X[i] - (Bx + Ax * Z[i]);
XY[i][2] = Y[i] - (By + Ay * Z[i]);
}
}
model {
sgm1 ~ uniform(0.001, 1000);
sgm2 ~ uniform(0.001,
1000);
rho ~ uniform(-0.999,
0.999);
Ax ~ normal(0.0, 1000.0);
Ay ~ normal(0.0, 1000.0);
Bx ~ normal(0.0, 1000.0);
By ~ normal(0.0, 1000.0);
for (i
in 1:N) {
XY[i] ~ multi_normal([0.0, 0.0],
[[sgm1*sgm1, rho*sgm1*sgm2],
[rho*sgm2*sgm1, sgm2*sgm2]]);
}
}
リスト２ リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプト（）
import csv
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import seaborn as sb
import scipy.stats
as ss
import arviz
as az
def calc_b_a(y,
x):
"""
単回帰分析
"""
X = np.vstack((np.ones(len(y)), x)).T
ba
= np.linalg.inv(X.T@X)@X.T@np.reshape(y,(len(y),1))
return ba[0][0],
ba[1][0]
def calc_r(x,
y):
"""
相関係数の計算
"""
num = np.sum((x
- np.mean(x)) * (y - np.mean(y)))
den = np.sum((x
- np.mean(x))**2) * np.sum((y
- np.mean(y))**2)
return num / (den**0.5)
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}
model = CmdStanModel(stan_file='pcor.stan')
fit = model.sample(data=Data)
print(fit.summary())
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata, var_names=['rho', 'sgm1', 'sgm2'])
plt.tight_layout()
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
plt.figure(figsize=(8,6))
plt.subplot(221)
sb.kdeplot(fit['Ax'])
plt.title('Posterior distribution of Ax')
plt.subplot(222)
sb.kdeplot(fit['Bx'])
plt.title('Posterior distiribution
of Bx')
plt.subplot(223)
sb.kdeplot(fit['Ay'])
plt.title('Posterior distribution of Ay')
plt.subplot(224)
sb.kdeplot(fit['By'])
plt.title('Posterior distribution of By')
plt.tight_layout()
plt.show()
sb.kdeplot(fit['rho'])
plt.xlabel('Partial Correlation', fontsize=12)
p_pos = np.mean(fit['rho']>0)
q1, q2, q3 = np.percentile(fit['rho'],
[25, 50, 75])
plt.title('Posterior distribution of Partial
Correlation' +
f'\nP(partial corr>0)={p_pos:.5f},'
+
f'\nQ1={q1:.3f}, Med.={q2:.3f}, Q3={q3:.3f}',
fontsize=14)
plt.tight_layout()
plt.show()
XZ = []
YZ = []
for i in
range(len(X)):
XZ.append(np.median(fit[f'XY[{i+1},1]']))
YZ.append(np.median(fit[f'XY[{i+1},2]']))
plt.figure(figsize=(9,8))
plt.subplot(221)
plt.scatter(Z,X)
plt.xlabel('Z', fontsize=14)
plt.ylabel('X', fontsize=14)
min_Z = np.min(Z)
max_Z = np.max(Z)
Bx_hat = np.median(fit['Bx'])
Ax_hat = np.median(fit['Ax'])
plt.plot([min_Z, max_Z],
[Bx_hat+Ax_hat*min_Z,
Bx_hat+Ax_hat*max_Z],
c='r')
r = calc_r(Z,X)
plt.title(f'r={r:.3f}',
fontsize=14)
plt.subplot(223)
plt.scatter(X,Y)
plt.xlabel('X', fontsize=14)
plt.ylabel('Y', fontsize=14)
b, a = calc_b_a(Y,X)
min_X = np.min(X)
max_X = np.max(X)
plt.plot([min_X, max_X],
[b+a*min_X, b+a*max_X],
c='r')
r = calc_r(Y,X)
plt.title(f'r={r:.3f}',
fontsize=14)
plt.subplot(224)
plt.scatter(Z,Y)
plt.xlabel('Z', fontsize=14)
plt.ylabel('Y', fontsize=14)
By_hat = np.median(fit['By'])
Ay_hat = np.median(fit['Ay'])
plt.plot([min_Z, max_Z],
[By_hat+Ay_hat*min_Z,
By_hat+Ay_hat*max_Z],
c='r')
r = calc_r(Z,Y)
plt.title(f'r={r:.3f}',
fontsize=14)
plt.subplot(222)
plt.scatter(XZ, YZ)
plt.xlabel('X.Z', fontsize=14)
plt.ylabel('Y.Z', fontsize=14)
b, a = calc_b_a(YZ,
XZ)
min_XZ = np.min(XZ)
max_XZ = np.max(XZ)
plt.plot([min_XZ, max_XZ],
[b+a*min_XZ, b+a*max_XZ], c='r')
r = np.median(fit['rho'])
plt.title(f'r={r:.3f}',
fontsize=14)
plt.tight_layout()
plt.show()
fout.close()
print('Results.txt was saved.')
Mardia, K.V., Kent, J.T., & Bibby, J.M. (1979) Multivariate analysis. Academic Press.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
岡本安晴（1998） Delphiで学ぶデータ分析法．CQ出版社