拙著「Pythonでデータ分析」の第６章「重回帰分析」の更新と補足
最小2乗法とベイズ分析
重回帰分析のスクリプト例を、最小2乗法に基づくものとベイズ分析を行うものとを作成してみた。
データをExcelファイルとして用意するスクリプトファイルをリスト１のように用意した。スクリプトファイルなどは、圧縮ファイルc06files.zipにまとめたので、ダウンロード展開して自由に利用されたい。
入力データは、Excelファイルとして図１のように用意する。
図１
第1行目に変数名を設定する。データは、２行目以降に書く。
第１列目には、ケースID用の値（文字列）を設定する。
第２列目に従属変数の値を設定する。
第３列目以降に、独立変数の値を設定する。
リスト１のスクリプトを実行すると、入力ファイル名を聞いてくる。
入力データファイル(*.xlsx) =
入力データファイル名（Excelであるので「*.xlsx」の形式である）、図１の場合は「予想最高気温.xlsx」を設定する。
入力データファイル(*.xlsx) = 予想最高気温.xlsx
入力データファイル名を設定すると、計算が始まり、（予測値、データ値）の散布図が表示される（図２）。
図２
図２のフォームを閉じると、図３のフォームが表示される。これは、変数（予測値、データ値）を標準化したときの散布図である（グラフの目盛りに注意）。
図３
図３のフォームを閉じると、スクリプトの実行終了となり、次のメッセージが表示されている。
Results.txt
was saved.
ファイルResults.txtを開くと、以下のような内容である。
入力データファイル = 予想最高気温.xlsx
切片
20.750
偏回帰係数
緯度
-0.872
西経
0.098
重相関係数 = 0.8146660647977868
多重決定係数 = 0.6636807971331118
標準化回帰係数
緯度
-0.815
(SD: 5.784)
西経
0.273 (SD:
17.235)
上記の出力において、切片は独立変数を中心化（centering）したときの値であることに注意（Gelmanら、2021、p.137）。すなわち、独立変数の平均値における回帰モデルの値である。北緯０°、西経０°における回帰モデルの値ではない。
独立変数は中心化のみを行っているので、偏回帰係数は緯度１°の差に対応する回帰モデルの値の差が−0.872℃であること、西経１°の差に対応する回帰モデルの値の差が0.098℃であることが示されている。
標準化回帰係数では、すべての値が標準化されているので、独立変数の値が１標準偏差異なるときの回帰モデルの値の標準偏差を単位としたときの差である。西経の偏回帰係数の値0.098が標準化回帰係数では0.273と大きくなっているが、これは独立変数「西経」の標準偏差が17.235と大きいこと、すなわち１標準偏差の表す西経の度の値が大きいことに対応するものである。
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
import pandas as pd
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import japanize_matplotlib
"""
Excelファイルからの読み込み
"""
in_file_name = input('入力データファイル(*.xlsx) = ')
xlsx = pd.ExcelFile(in_file_name)
data_xlsx = pd.read_excel(xlsx)
"""
出力ファイルの設定
"""
f = open('Results.txt', 'w')
f.write('\n入力データファイル = ' + in_file_name)
keys = data_xlsx.keys()
# １行目の変数名の読み込み
print('Keys...\n', keys)
case_name = keys[0]
case_id = [ v for v in
data_xlsx[case_name]]
print('ケースID:\n', case_id)
y_name = keys[1]
y = np.array([v for v in
data_xlsx[y_name]]) # 従属変数yの読み込み
print('y =\n', y)
print('Mean of y =', np.mean(y))
print('SD of y =', np.std(y))
y
N = len(y)
# データ数
print('N =', N)
M = len(keys) - 2 # 独立変数の数
print('M =', M)
keys = np.array(keys)
x_names = keys[2:] # 独立変数名リスト
print('x_names =', x_names)
X = []
# 独立変数
for x_name in x_names:
X.append(data_xlsx[x_name].values)
X = np.array(X)
print('X =\n', X)
Means_X = np.mean(X, axis=1)
print('Means =\n', Means_X)
SDs = np.std(X, axis=1)
print('SDs =\n', SDs)
Xc = X.T - np.mean(X, axis=1) # 独立変数の中心化
Zx = Xc / SDs
# 独立変数の標準化
"""
Xc <- [1, Xc] 参照（6.1.3）式
"""
Xc = np.concatenate((np.ones((N,1)),
Xc), axis=1)
print('Xc =\n', Xc)
b = np.linalg.inv(Xc.T @ Xc) @ Xc.T @
y # （6.1.5）式
print(f'{"切片":>20s}{b[0]:>20.3f}')
print('偏回帰係数')
for i in range(M):
print(f'{x_names[i]:>20s}{b[i+1]:>20.3f}')
f.write(f'\n\n{"切片":>20s}{b[0]:>20.3f}')
f.write('\n\n偏回帰係数')
for i in range(M):
f.write(f'\n{x_names[i]:>20s}{b[i+1]:>20.3f}')
y_hat = Xc @ b
# 回帰モデルによる予測値
R = ss.pearsonr(y_hat, y)[0] # 重相関係数
print('\n重相関係数 =', R)
f.write(f'\n\n重相関係数 = {R}')
R2 = R**2
# 多重決定係数
print('多重決定係数 =', R2)
f.write(f'\n多重決定係数 = {R2}')
"""
散布図（y_hat, y）
"""
plt.plot(y_hat, y, 'o')
for i in range(N):
plt.text(y_hat[i], y[i],
case_id[i])
plt.plot([np.min(y_hat),
np.max(y_hat)], [np.min(y_hat), np.max(y_hat)],
label = 'y = $\hat{y}$')
plt.xlabel('$\hat{y}$', fontsize = 16)
plt.ylabel('y', fontsize = 16)
plt.title('Scatter Plot ($\hat{y},
y$)', fontsize = 18)
plt.legend(fontsize = 14)
plt.savefig('FigScatter.png')
plt.show()
"""
標準化得点による重回帰分析
"""
y_mean = np.mean(y)
y_sd = np.std(y)
yz = (y - y_mean) / y_sd
print('\ny_mean =', y_mean)
print('y_sd =', y_sd)
print('yz =\n', yz)
bz = np.linalg.inv(Zx.T @ Zx) @ Zx.T @
yz # 拙著 p.109参照
print('\n標準化回帰係数\n')
for i in range(M):
print(f'{x_names[i]:>20s}{bz[i]:>20.3f}{"(SD:
":>20s}{SDs[i]:.3f})')
f.write(f'\n\n標準化回帰係数\n')
for i in range(M):
f.write(f'\n{x_names[i]:>20s}{bz[i]:>20.3f}{"(SD:
":>20s}{SDs[i]:.3f})')
yz_hat = Zx @ bz
"""
散布図（y_hat, y）
"""
plt.plot(yz_hat, yz, 'o')
for i in range(N):
plt.text(yz_hat[i], yz[i],
case_id[i])
plt.plot([np.min(yz_hat),
np.max(yz_hat)], [np.min(yz_hat), np.max(yz_hat)],
label = 'y = $\hat{y}$')
plt.xlabel('$\hat{y}$', fontsize = 16)
plt.ylabel('y', fontsize = 16)
plt.title('Scatter Plot ($\hat{y},
y$)', fontsize = 18)
plt.legend(fontsize = 14)
plt.savefig('FigScatterZ.png')
plt.show()
f.close()
print('\nResults.txt was saved.')
重回帰モデルのベイズ分析例は、このウェブサイトで説明している。