偏相関係数
疑似相関
CmdStanPy対応：2025.11、2026.02
2つの変数XとYの相関係数が、第３の変数Zの影響を除くと相関が無くなることがある。この疑似相関（spurious correlation）に対して、第３の変数Zの影響を単回帰モデルによって変数XとYから除き、その残りの部分について相関係数を求めたものは偏相関係数（partial correlation coefficient）と呼ばれている（Mardia et al., 1979）。すなわち、２つの単回帰モデル
において、
と
の相関係数を偏相関係数と呼び、
で表す。
である。
XとY、YとZ、XとZの間の相関係数を
、
、
で表すと次式が成り立つ。
＜岡本安晴「Delphiで学ぶデータ分析法」CQ出版社＞のpp．32-35に証明がある。
式（４）を導く考え方については、＜岡本安晴「データを使いこなすための統計入門」22世紀アート、電子書籍＞のp.40-41でも説明している。
ここでは、式（４）により代数的に偏相関係数を求めるのではなく、ベイズ分析により偏相関係数の事後分布を求める。事後分布から、偏相関係数がどの範囲にあるかの推測を行うことができる。
定義より、偏相関係数を
と置くとき、式（１）と式（２）に対して
と置くことができる。
上式に基づき、リスト1に示すStanスクリプトを用意した。
リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプトをリスト3に用意した。
ファイルは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スクリプトファイル（est_cors.py）、および入力データファイ（DataXYZ.csv）を同じフォルダに置く。リスト３のest_cors.pyは、3変数X、Y、Zの相関行列のベイズ分析を行うStanスクリプトest_r.stan（リスト２）も用いているので、est_t.stanファイルも同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境で、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイト、あるいはこのウェブサイトなどに用意した。
(stan) ****/SpurCorrFiles$ python est_cors.py
上のコマンドを実行すると、入力データファイル名の設定が求められる。
(stan) ****/SpurCorrFiles$
python est_cors.py
Input file name (*.csv) = DataXYZ.csv
上の例では、図１のファイル名が設定されている。
ファイル名が設定されると、まず、リスト２のStanスクリプトにより、３変数間の相関行列のベイズ分析が行われる。
相関行列は、コレスキー分解（岡本、2009、pp.46-47）を行い、コレスキー因子に対してサンプリングを行っている。
Rs = L * (L');
L ~ lkj_corr_cholesky(1.0);
リスト２のStanスクリプトによるMCMCサンプリングのトレース図は、図A.1のように表示される。
図A.1
図A.1のWindowを閉じると、相関係数の結果が図A.2のように表示される。
図A.2
図A.2のWindowを閉じると、リスト１のStanスクリプトによる偏相関係数のベイズ分析が始まる。
次のコードで回帰モデルが設定されている。
XY[i][1] = X[i] - (Bx + Ax * Z[i]);
XY[i][2] = Y[i] - (By + Ay * Z[i]);
回帰モデルの残差に対して、相関係数が次のコードで求められている。
XY[i] ~ multi_normal([0.0, 0.0],
[[sgm1*sgm1, rho*sgm1*sgm2],
[rho*sgm2*sgm1, sgm2*sgm2]]);
MCMCサンプリングが終わると、トレース図が表示される（図A.3）。
図A.3
図A.３のWindowを閉じると、図A.４のグラフが表示される。
図A.４
回帰直線のパラメータの事後分布のグラフである。
図A４のWindowを閉じると、図A.５のグラフが表示される。
図A.５
偏相関係数の事後分布のグラフである。偏相関係数が正である事後確率は、0.31375であると表示されている。事後分布は、ほぼ原点を中心にして分布していると言える。
図A.５のWindowを閉じると、図A.６のグラフが表示される。
図A.６
点（Z,X）の散布図が左上、点（X.Z、Y.Z）の散布図が右上に示されている。X.ZとY.Zは、式（１）における残差と式（２）における残差を表している。下段の散布図は、点（X,Y）の散布図が左側、点（Z,Y）の散布図が右側に示されている。散布図における赤の実線は、回帰直線である。散布図の上に表示されているｒの値は相関係数である。
変量XとYの相関係数はr=0.496であるが、偏相関係数はr=0.035であり、図A.５の事後分布から、偏相関係数は０であると言えそうである。
図A.６のWindowを閉じると、スクリプトの実行終了である。
次に、図２のデータを分析してみる。
リスト１、リスト２、リスト３のスクリプトファイル、および図２のデータファイルを同じフォルダに置いて、そのフォルダにカレントディレクトリを移し、CmdStanPyのインストールされている環境において、次のコマンドを実行する。
(stan) ****/SpurCorrFiles$ python est_cors.py
入力データファイル名を聞いてくるので、図２のデータファイル名を設定する。
(stan) ****/SpurCorrFiles$
python est_cors.py
Input file name (*.csv) = DataXYZW.csv
まず、３変数間の相関行列のベイズ分析が始まる。
MCMCサンプリングのトレース図が、図B.1のように表示される。
図B.1
図B.1のWindowを閉じると、図B.2のグラフが表示される（図B.2）。
図B.2
図B.2のWindowを閉じると、リスト１のStanスクリプトによる偏相関係数のベイズ分析が始まる。
MCMCサンプリングが終了すると、トレース図が表示される（図B.3）。
図B.３
図B.3のWindowを閉じると、図B.4のグラフが表示される。
図B.4
図B.4のWindowを閉じると、図B.5のグラフが表示される。
図B.5
偏相関係数の事後分布のグラフであるが、偏相関係数が正である事後確率は、0.00000と表示されている。事後分布の中央値は、―0.514である。偏相関係数は、負であると言えそうである。
図B.5のWindowを閉じると、図B.6のグラフが表示される。
図B.6
左下の（X,Y）の散布図では、相関係数は0.629であり、回帰直線は右上がりの直線である。変量Zの影響を除いた点（X.Z、Y.Z）の散布図では、偏相関係数はー0.514で、回帰直線は右下がりの直線である。
図B.6の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;
/*
Zの影響を単回帰モデルにより除いている。
*/
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]]);
}
}
リスト２ 3変数X、Y、Zの相関行列を求めるStanスクリプト（est_r.stan）
data {
int K;
int n;
array[n] real x;
array[n] real y;
array[n] real z;
}
parameters {
real mu1;
real mu2;
real mu3;
vector[K] SDs;
cholesky_factor_corr[K]
L;
}
transformed parameters {
real v;
matrix[K,K] Rs;
matrix[K,K] Sigma;
/*
Cholesky factor L から相関係数行列Rを算出
R = L(L')
*/
Rs = L * (L');
/*
相関係数rijと標準偏差sdiおよびsdjから共分散Sigmaijを算出
Sigmaij = sdi * rij * sdj
*/
for (i
in 1:K) {
for
(j in 1:i) {
Sigma[i,j] = SDs[i]
* Rs[i,j] * SDs[j];
Sigma[j,i] = Sigma[i,j];
}
}
}
model {
L ~ lkj_corr_cholesky(1.0);
SDs[1] ~ uniform(0.0001,
1000);
SDs[2] ~ uniform(0.0001,
1000);
SDs[3] ~ uniform(0.0001,
1000);
for (i
in 1:n) {
[x[i],y[i],z[i]]
~ multi_normal([mu1,mu2,mu3], Sigma);
}
}
リスト３ リスト１のStanスクリプトを用いて偏相関係数を求めるPythonスクリプト（est_cors.py）
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
import networkx
as nx
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 = {'K':3, 'n': len(X), 'x':X, 'y':Y, 'z':Z}
#
#
相関行列のベイズ推定
#
sm = CmdStanModel(stan_file='est_r.stan')
fit = sm.sample(data=Data,
iter_warmup=10_000)
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data, var_names=['mu1', 'mu2', 'mu3', 'SDs', 'Rs'],
figsize=(12,7))
# トレースプロット
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
print('Keys...\n', dframe.keys())
plt.figure(figsize=(12,7))
plt.subplot(221)
sb.kdeplot(dframe['Rs[1,2]'])
plt.xticks([-1,-0.5,0,0.5,1])
plt.xlabel('$r_{XY}$', fontsize=14)
plt.subplot(223)
sb.kdeplot(dframe['Rs[1,3]'])
plt.xticks([-1,-0.5,0,0.5,1])
plt.xlabel('$r_{XZ}$', fontsize=14)
plt.subplot(224)
sb.kdeplot(dframe['Rs[2,3]'])
plt.xticks([-1,-0.5,0,0.5,1])
plt.xlabel('$r_{YZ}$', fontsize=14)
plt.subplot(222)
G = nx.Graph()
G.add_edge('X','Y')
G.add_edge('X','Z')
G.add_edge('Y','Z')
med12 = np.median(dframe['Rs[1,2]'])
med13 = np.median(dframe['Rs[1,3]'])
med23 = np.median(dframe['Rs[2,3]'])
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels
= True, font_size=20, # グラフの描画
node_size=1200, node_color='orange',
edgecolors='m')
nx.draw_networkx_edge_labels(G, pos,
edge_labels={('X','Y'):f"r={med12:.3f}",
('X','Z'):f"r={med13:.3f}",
('Y','Z'):f"r={med23:.3f}"},
font_size=20)
plt.tight_layout()
plt.show()
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でデータ分析．丸善出版
岡本安晴（2009）統計学を学ぶための数学入門［下］、培風館
岡本安晴（1998） Delphiで学ぶデータ分析法．CQ出版社