２章 複数の変数の間の関係を見る
重回帰分析
物理の成績と数学の成績および体重との関係を調べるために図１のようなデータを収集したとする。
図１ Data2_1.xlsx
図１のデータは＜岡本安晴「心理学データ分析と測定」勁草書房＞の表2.1.1（仮想データ）から２０人分を抜き出したものである。１行目に変数名が書かれている。１列目はデータの識別値（例えば通し番号、数値でも文字列でも何でもよい）が設定され、２列目が従属変数、３列目からが独立変数である。
いま、物理、数学、体重の間に次式（2.1.1）の関係を想定する。式番号なども上記＜岡本安晴「心理学データ分析と測定」勁草書房＞におけるものを用いる。
まず、図１のデータを読み込むために、次のスクリプトを実行する。
flnm = input('Data File Name (*.xlsx) = ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
PandasのDataFrame型に読み込まれたデータを表すdataから、以下のスクリプトで独立変数値と従属変数値を取り出しXとYで表す。ここで、Xの第１列には１が設定されていることに注意。
VarLabels = ['Const']
+ [v for v in data.columns[2:]] # 独立変数名
print(VarLabels)
p = len(VarLabels)
# 独立変数の数
print('p = ', p)
n = len(data.values)
# データ数
print('n = ', n)
YLabel = data.columns[1]
# 従属変数名
Y = data[YLabel].values.reshape((n,1)) # 従属変数値
print(Y)
X = data.values[:,
2:]
X = np.hstack((np.array([1]*n).reshape((n,1)), X)) # 独立変数値
偏回帰係数
は次式
で求められるが、これは次のコードで算出される。
b = np.linalg.inv(X.T @ X) @ X.T @ Y
# 偏回帰係数の計算
信頼区間の計算では、まず次の行列
（p．25）が計算される。
この行列の対角成分の平方根
が求められ、
信頼区間が次式で与えられる。
スクリプトは以下のようである。ただし、
である。
Yhat = X @ b
SSE = ((Yhat
- Y) ** 2).sum()
MSE = SSE / (n - p)
s2_b = MSE * np.linalg.inv(X.T
@ X)
s_b = np.diag(s2_b)
** 0.5
t975 = ss.t.ppf(0.975,
n - p)
print()
for i in
range(p):
print(VarLabels[i] + ':')
print(' b[{0:}] = {1:.3f} 95%CI = [{2:.3}, {3:.3}]'.format(i, b[i][0],
b[i][0] - t975 * s_b[i],
b[i][0] + t975 * s_b[i]))
多重決定係数
と重相関係数
は、次式で算出される。
ここで、
である。これより、
である。
スクリプトは、以下のようである。
Ymean = Y.mean()
print('\nY<{0:}>_mean
= {1:}'.format(YLabel, Ymean))
SSReg = ((Yhat - Ymean) ** 2).sum()
SSTot = ((Y - Ymean)
** 2).sum()
R2 = SSReg /
SSTot
print('\nR^2 = ', R2)
R = R2 ** 0.5
print('R = ', R)
完全なプログラムをリスト１に示す。スクリプトファイルなどは圧縮ファイルfiles2.zipとしてまとめたので、ダウンロードすれば自由に利用できる。
リスト１のスクリプトを実行すると、入力データファイル名の設定が求められる（図２）。
図２
例えば、図１のデータファイルの名前を設定すると計算が始まり、計算結果が図３のように表示される。
図３
独立変数を標準得点化する場合
独立変数を標準得点に変換する場合は、次のスクリプトで簡単にできる。
X = ss.zscore(X)
# 独立変数の標準得点化
この場合のPythonスクリプトをリスト２に示した。このスクリプトファイルも圧縮ファイルfiles2.zipにまとめたので、ダウンロードして自由に利用できる。
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import scipy.stats
as ss
flnm = input('Data File Name (*.xlsx) = ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
VarLabels = ['Const']
+ [v for v in data.columns[2:]] # 独立変数名
print(VarLabels)
p = len(VarLabels)
# 独立変数の数
print('p = ', p)
n = len(data.values)
# データ数
print('n = ', n)
YLabel = data.columns[1]
# 従属変数名
Y = data[YLabel].values.reshape((n,1)) # 従属変数値
print(Y)
X = data.values[:,
2:]
X = np.hstack((np.array([1]*n).reshape((n,1)), X)) # 独立変数値
print(X)
b = np.linalg.inv(X.T
@ X) @ X.T @ Y
# 偏回帰係数の計算
#
#
95%信頼区間の計算
#
Yhat = X @ b
SSE = ((Yhat
- Y) ** 2).sum()
MSE = SSE / (n - p)
s2_b = MSE * np.linalg.inv(X.T
@ X)
s_b = np.diag(s2_b)
** 0.5
t975 = ss.t.ppf(0.975,
n - p)
print()
for i in
range(p):
print(VarLabels[i] + ':')
print(' b[{0:}] = {1:.3f} 95%CI = [{2:.3}, {3:.3}]'.format(i, b[i][0],
b[i][0] - t975 * s_b[i],
b[i][0] + t975 * s_b[i]))
#
# 多重決定係数と重相関係数の計算
#
Ymean = Y.mean()
print('\nY<{0:}>_mean
= {1:}'.format(YLabel, Ymean))
SSReg = ((Yhat - Ymean) ** 2).sum()
SSTot = ((Y - Ymean)
** 2).sum()
R2 = SSReg /
SSTot
print('\nR^2 = ', R2)
R = R2 ** 0.5
print('R = ', R)
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import scipy.stats
as ss
flnm = input('Data File Name (*.xlsx) = ')
xlsx = pd.ExcelFile(flnm)
# 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
VarLabels = ['Const']
+ [v for v in data.columns[2:]] # 独立変数名
print(VarLabels)
p = len(VarLabels)
# 独立変数の数
print('p = ', p)
n = len(data.values)
# データ数
print('n = ', n)
YLabel = data.columns[1] # 従属変数名
Y = data[YLabel].values.reshape((n,1)) # 従属変数値
print(Y)
X = data.values[:,
2:]
print(X)
X = ss.zscore(X)
# 独立変数の標準得点化
print(X)
X = np.hstack((np.array([1]*n).reshape((n,1)), X)) # 独立変数値
print(X)
b = np.linalg.inv(X.T
@ X) @ X.T @ Y
# 偏回帰係数の計算
#
#
95%信頼区間の計算
#
Yhat = X @ b
SSE = ((Yhat
- Y) ** 2).sum()
MSE = SSE / (n - p)
s2_b = MSE * np.linalg.inv(X.T
@ X)
s_b = np.diag(s2_b)
** 0.5
t975 = ss.t.ppf(0.975,
n - p)
print()
for i in
range(p):
print(VarLabels[i] + ':')
print(' b[{0:}] = {1:.3f} 95%CI = [{2:.3}, {3:.3}]'.format(i, b[i][0],
b[i][0] - t975 * s_b[i],
b[i][0] + t975 * s_b[i]))
#
# 多重決定係数と重相関係数の計算
#
Ymean = Y.mean()
print('\nY<{0:}>_mean
= {1:}'.format(YLabel, Ymean))
SSReg = ((Yhat - Ymean) ** 2).sum()
SSTot = ((Y - Ymean)
** 2).sum()
R2 = SSReg /
SSTot
print('\nR^2 = ', R2)
R = R2 ** 0.5
print('R = ', R)