重回帰分析
独立変数が
と
の２個である仮想データ
例を図１に示す。
図１
図１に対する重回帰モデルを以下のように設定する。
パラメータの事前分布は、平らな事前分布として標準偏差の大きい正規分布と指数分布を以下のように設定する。
上記の確率モデルによるベイズ分析を行うためにStanスクリプトをリスト１のように用意した。
スクリプトファイルとデータファイルは、圧縮ファイルfiles_mregress.zipにまとめた。ダウンロード・展開すれば自由に利用できる。
Stanについては、岡本（2019）で説明している。
リスト１ 重回帰分析のStanスクリプト（mv_regress.stan）
data {
int
N;
real Y[N];
real X1[N];
real X2[N];
}
parameters {
real b0;
real b1;
real b2;
real<lower = 0.0> sgm;
}
model {
b0 ~ normal(0, 1000);
b1 ~ normal(0, 1000);
b2 ~ normal(0, 1000);
sgm
~ exponential(0.01);
for (i
in 1:N)
Y[i] ~ normal(b0 + b1 * X1[i] + b2
* X2[i], sgm);
}
リスト１のStanスクリプトを使って図１のデータの重回帰分析を行うPythonスクリプトをリスト２のように用意した。
リスト２ 重回帰分析を行うPythonスクリプト（MRegress.py）
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import pystan
from matplotlib
import cm
import seaborn
as sb
df = pd.read_csv('MVData.csv')
data = df.values
print('data =')
Y = data.T[0]
X1 = data.T[1]
X2 = data.T[2]
for y, x1, x2 in zip(Y, X1, X2):
print(y, x1, x2)
sm = pystan.StanModel(file
= 'mv_regress.stan')
print('Compiling ende.')
Data = {'N':len(Y),
'Y':Y, 'X1':X1, 'X2':X2}
fit = sm.sampling(data
= Data, n_jobs = 1)
print(fit)
b0_hat = np.median(fit['b0'])
b1_hat = np.median(fit['b1'])
b2_hat = np.median(fit['b2'])
sb.kdeplot(fit['b0'])
plt.title('b0: median =
{0:.2f}'.format(b0_hat))
plt.savefig('Figb0.png')
plt.show()
sb.kdeplot(fit['b1'])
plt.title('b1: median = {0:.2f}'.format(b1_hat))
plt.savefig('Figb1.png')
plt.show()
sb.kdeplot(fit['b2'])
plt.title('b2: median =
{0:.2f}'.format(b2_hat))
plt.savefig('Figb2.png')
plt.show()
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(X1, X2, Y, color='blue')
Xs = np.linspace(min(X1),
max(X1), 100)
Ys = np.linspace(min(X2),
max(X2), 100)
Xg, Yg = np.meshgrid(Xs, Ys)
Z = np.empty((100,
100))
for iy in
range(100):
for ix in range(100):
Z[iy][ix] = b0_hat + b1_hat * Xs[ix]
+ b2_hat * Ys[iy]
ax.plot_surface(Xg, Yg, Z, cmap = cm.viridis,
alpha = 0.3)
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_zlabel('Y')
ax.set_title('Multivariate Regression')
for i in
range(len(Y)):
plt.plot([X1[i],X1[i]], [X2[i],X2[i]],
[b0_hat + b1_hat*X1[i] + b2_hat*X2[i], Y[i]], c = 'm')
plt.plot([0], [0], [b0_hat], 'D', color = 'r',
label = (r'(0,0,$\hat{b}_0}$)') + '\n' +
(r'$\hatb_0$={0:.2f}'.format(b0_hat)))
plt.legend()
plt.title('You can rotate this graph by
dragging')
plt.savefig('FigRegress.png')
plt.show()
リスト２を実行すると、以下の結果を得る。
図２はb0の事後分布である。
図２ b0の事後分布
図３は、b1の事後分布である。
図３ b1の事後分布
図４は、b2の事後分布である。
図４ b2の事後分布
図５は、回帰平面である。
図５ 回帰平面
図５では、重回帰モデルの定数項b0は点
における値であるので、データ点の集まりから離れた位置に
がプロットされている。
独立変数を中心化すると、定数項は独立変数の平均値での重回帰平面の値となる。
偏回帰係数b1、b2の値は、それぞれの独立変数の値が１だけ異なるときの回帰平面の値の違いを表す。これを独立変数の値が標準偏差だけ異なるときの回帰平面の値の違いを表すようにするときは、それぞれの独立変数の値をその標準偏差で割ったものを用いればよい。
変数
の中心化と標準偏差を１に揃えることは変量の標準化と呼ばれているが、標準化後の変量を
で表すと、次式で与えられる。
ここで、
は平均値、
は標準偏差を表す。
独立変数の標準化を行って、重回帰分析を行うPythonスクリプトをリスト３に示す。
リスト３ 独立変数の標準化による重回帰分析（MRegress_Z.py）
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import pystan
import seaborn
as sb
from matplotlib
import cm
df = pd.read_csv('MVData.csv')
data = df.values
Y = data.T[0]
X1 = data.T[1]
X2 = data.T[2]
for y, x1, x2 in zip(Y, X1, X2):
print(y, x1, x2)
Z1 = (X1 - np.mean(X1))
/ np.std(X1)
Z2 = (X2 - np.mean(X2))
/ np.std(X2)
print('X1: mean = {0:.3f} std =
{1:.3f}'.format(np.mean(X1), np.std(X1)))
print('X2: mean = {0:.3f} std =
{1:.3f}'.format(np.mean(X2), np.std(X2)))
sm = pystan.StanModel(file
= 'mv_regress.stan')
print('Compiling ende.')
Data = {'N':len(Y),
'Y':Y, 'X1':Z1, 'X2':Z2}
fit = sm.sampling(data
= Data, n_jobs = 1)
print(fit)
b0_hat = np.median(fit['b0'])
b1_hat = np.median(fit['b1'])
b2_hat = np.median(fit['b2'])
sb.kdeplot(fit['b0'])
plt.title('b0: median =
{0:.2f}'.format(b0_hat))
plt.savefig('Figb0_c.png')
plt.show()
sb.kdeplot(fit['b1'])
plt.title('b1: median =
{0:.2f}'.format(b1_hat))
plt.savefig('Figb1_c.png')
plt.show()
sb.kdeplot(fit['b2'])
plt.title('b2: median = {0:.2f}'.format(b2_hat))
plt.savefig('Figb2_c.png')
plt.show()
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(Z1, Z2, Y, color='blue')
Xs = np.linspace(min(Z1),
max(Z1), 100)
Ys = np.linspace(min(Z2),
max(Z2), 100)
Xg, Yg = np.meshgrid(Xs, Ys)
Z = np.empty((100,
100))
for iy in
range(100):
for ix in range(100):
Z[iy][ix] = b0_hat + b1_hat * Xs[ix]
+ b2_hat * Ys[iy]
ax.plot_surface(Xg, Yg, Z, cmap = cm.viridis,
alpha = 0.3)
ax.set_xlabel('Z1')
ax.set_ylabel('Z2')
ax.set_zlabel('Y')
ax.set_title('Multivariate Regression')
for i in
range(len(Y)):
plt.plot([Z1[i],Z1[i]], [Z2[i],Z2[i]],
[b0_hat + b1_hat*Z1[i] + b2_hat*Z2[i], Y[i]], c = 'm')
plt.plot([0], [0], [b0_hat], 'D', color = 'r',
label = (r'(0,0,$\hat{b}_0}$)') + '\n' +
(r'$\hatb_0$={0:.2f}'.format(b0_hat)))
plt.legend()
plt.title('You can rotate this graph by
dragging')
plt.savefig('FigRegress_c.png')
plt.show()
リスト３を実行すると、以下の結果を得る。
定数項b0の事後分布は、図６のようである。
図６ b0の事後分布
偏回帰係数b1の事後分布は、図７のようになる。
図７ b1の事後分布
偏回帰係数b2の事後分布は、図８のようである。
図8 b2の事後分布
重回帰平面は、図９のようである。
図９ 回帰平面
定数項の表す点
が、データの独立変数Z1,Z2の平均値０に対応する位置の重回帰平面の値であるので、データ全体の中央に位置している。独立変数Z1、Z2の変域も標準化されているので同じ範囲であるといえる。