重回帰分析
重回帰分析（multiple linear regression analysis）のプログラムを＜岡本安晴「心理学データ分析と測定」勁草書房＞の第2章の式に基づいて作成した。プログラムファイルとサンプルデータファイルは、圧縮ファイルMRgrsnPrg.zipにまとめた。
重回帰モデルは
と表せる。式（１）のモデルを、例えば表2.1.1のデータに当てはめる。
|
表2.1.1 （物理、数学、身長）の50人分のデータ |
|
(53, 73, 167) (65, 53, 152) (54, 50, 145) (59, 58, 153) (58, 49, 176) (51, 52, 148) (52, 49, 159) (63, 69, 153) (64, 63, 158) (75, 83, 161) (66, 58, 171) (59, 57, 154) (75, 90, 172) (65, 76, 167) (70, 61, 160) (62, 60, 152) (58, 44, 166) (79, 98, 160) (65, 65, 161) (63, 54, 161) (60, 47, 173) (64, 67, 165) (67, 71, 168) (52, 53, 165) (60, 70, 162) (63, 70, 171) (69, 72, 173) (57, 58, 159) (58, 65, 165) (51, 51, 157) (60, 63, 163) (50, 42, 151) (52, 43, 177) (67, 69, 149) (59, 67, 167) (44, 35, 168) (79, 85, 157) (54, 66, 180) (75, 82, 150) (60, 74, 156) (62, 63, 165) (54, 48, 164) (69, 73, 166) (74, 87, 163) (72, 81, 160) (56, 66, 165) (70, 77, 174) (64, 60, 166) (71, 72, 145) (70, 69, 151) |
プログラムでは、最小2乗法の基準で（１）式の係数を求めている。
プログラムを実行すると、入力データファイル名の設定が求められる（図１）。
図１
入力データファイルは、図２のような形式で用意する。
・
・
・
図２
第１文字がスラッシュ / で始まる行の次の行に、独立変数の数を書く。この第１文字がスラッシュで始まる行までは、ファイルは読み飛ばされるので、コメントなどを自由に書くことができる。
独立変数の数を書いた次の行に、従属変数名と独立変数名をタブ、あるいは半角空白文字で区切って書き並べる。
次の行から、１行に１組（従属変数と独立変数）のデータを、上の変数名の順序に合わせて書いていく。最後のデータの組を書いたら、その次の行は第１文字がスラッシュ / である行を置いて、データの終わりであることを示す。
データファイルが、プログラムファイルと同じフォルダ内にあるときは、データファイル名のみを入力すればよい（図３）。データファイル名を設定してEnterキーを押すと、次に出力ファイル名の設定が求められる。
図３
出力ファイル名はテキストファイル名（拡張子が「.txt」）を設定する（図４）。出力ファイル名の設定後、Enterキーを押すと計算が始まる。
図４
計算終了後、出力ファイルを開くと、表2.1.1のデータの場合、以下のようになっている。
Input data file =
Data211.txt
Labels...
Y: Physics
X[0]: (Const)
X[1]: Math
X[2]: Height
Input data...
53.0 73.0 167.0
65.0 53.0 152.0
54.0 50.0 145.0
59.0 58.0 153.0
58.0 49.0 176.0
51.0 52.0 148.0
・
・
・
72.0 81.0 160.0
56.0 66.0 165.0
70.0 77.0 174.0
64.0 60.0 166.0
71.0 72.0 145.0
70.0 69.0 151.0
b[0]<(Const)> =
43.0634604
b[1]<Math> =
0.4944112
b[2]<Height> =
-0.0766585
Physics = 43.0634604 +
0.4944112 * Math + -0.0766585 * Height + Residual
R2 = 0.6802908
R = 0.8247974
b[0]<(Const)>:
t =
3.2620618 df = 47 p-value = 0.0020633
b[1]<Math>:
t =
9.9499663 df = 47 p-value = 0.0000000
b[2]<Height>:
t =
-0.9704341 df = 47 p-value = 0.3367980
95% Confidence Intervals:
b[0]<(Const)>:
[16.50587, 69.62105]
b[1]<Math>:
[0.39445, 0.59437]
b[2]<Height>:
[-0.23557, 0.08226]
まず、入力データファイル名が出力されている。続いて、変数名が出力され、その後、読み込まれたデータ値が出力されている。これらにより、データが正しく読み込まれたかどうかのチェックができる。
読み込みデータの出力後、算出された
の推定値がb[j]の値として変数名とともに出力され、最後に重回帰モデル（１）の係数にb[j]の値を代入した形で式が書かれている。
式に続いて、多重決定係数
と重相関係数Rの値が、Ｒ2およびＲの値として書かれている。
次に、各係数の値のｔ値、自由度、およびｐ値（両側検定）が出力されている。
最後に、各係数の９５％信頼区間が書き出されている。
メインプログラムのリスト（スクリプト）を以下に示す。メインプログラムを含めたプログラムのソースコードファイル（スクリプト）とサンプルデータファイルは、圧縮ファイルMRgrsnPrg.zipとしてまとめた。
from math import *
import numpy as np
import numpy.linalg
as nl
from splitstring
import *
from calc_p_from_t
import *
from calc_t_from_p
import *
def CheckStr(str):
ck = False
if len(str) > 0:
if
str[0] == '/':
ck = True
return ck
#
# Prepare
the input data file
#
s_in =
input("Input data file = ")
f_in = open(s_in,
"r")
#
# Set the
contents of the input data file in the object data
#
data_f = f_in.readlines()
f_in.close()
#
# Prepare
the output file
#
s_out =
input("Output data file = ")
f_out = open(s_out,
"w")
f_out.write("Input
data file = " + s_in + "\n\n")
pos = 0
while not
CheckStr(data_f[pos]):
print(data_f[pos])
pos += 1
pos += 1
splitStr = CSplit('
', '\t', '\n')
splitStr.set_string(data_f[pos])
p =
int(splitStr.seek_nonempty(0))
print("Number of
independent variables, p = {}\n".format(p))
p += 1
pos += 1
splitStr.set_string(data_f[pos])
YLabel = splitStr.seek_nonempty(0)
XLabel = ['(Const)']
for i in range(1, p):
XLabel.append(splitStr.seek_nonempty(i))
f_out.write('\nLabels...\n')
f_out.write('Y: {}\n'.format(YLabel))
for i in range(p):
f_out.write('X[{0}]:
{1}\n'.format(i, XLabel[i]))
y = []
x = []
n = 0
pos += 1
while not
CheckStr(data_f[pos]):
splitStr.set_string(data_f[pos])
y.append(float(splitStr.seek_nonempty(0)))
x.append([1])
for j in range(1, p):
x[n].append(float(splitStr.seek_nonempty(j)))
n += 1
pos += 1
Y = np.matrix(y)
Y = np.transpose(Y)
X = np.matrix(x)
f_out.write("
")
f_out.write("Input
data...\n")
for i in range(n):
f_out.write("{0:10}".format(Y[i,0]))
for j in range(1, p):
f_out.write("{0:7}".format(X[i, j]))
f_out.write('\n')
XpXinv =
nl.inv(np.transpose(X) * X)
b = XpXinv *
np.transpose(X) * Y
f_out.write('\n')
for j in range(p):
f_out.write("b[{0}]<{1}> = {2:.7f}\n".format(j,
XLabel[j], b[j, 0]))
f_out.write('\n')
eq_str = YLabel +
" = " + ("{0:.7f}".format(b[0, 0]))
for j in range(1, p):
eq_str += " + " +
"{0:.7f} * {1}".format(b[j, 0], XLabel[j])
eq_str += " +
Residual"
f_out.write(eq_str +
'\n')
sumY = 0.0
for i in range(n):
sumY += Y[i, 0]
meanY = sumY / n
Yhat = X * b
SSReg = 0.0
SSTot = 0.0
for i in range(n):
SSReg += (Yhat[i, 0] -
meanY) ** 2.0
SSTot += (Y[i, 0] - meanY)
** 2.0
R2 = SSReg / SSTot
R = sqrt(R2)
f_out.write('\nR2 =
{0:.7}\n'.format(R2))
f_out.write('R = {0:.7}\n'.format(R))
SSE = np.transpose(Y
- X * b) * (Y - X * b)
MSE = SSE[0, 0] / (n
- p)
s2_b = MSE * XpXinv
t = []
for j in range(p):
t.append(b[j, 0] /
sqrt(s2_b[j,j]))
calc_p = C_cum_t(n -
p)
p_v = []
for j in range(p):
tp = calc_p.v(t[j])
if tp < 0.5:
p_v.append(2.0 * tp)
else:
p_v.append(2.0 * (1.0 - tp))
f_out.write('\n')
for j in range(p):
f_out.write('b[{0}]<{1}>:\n t = {2:.7f} df = {3} p-value = {4:.7f}\n'.\
format(j, XLabel[j], t[j], n - p, p_v[j]))
calc_t = CCalc_t(n - p)
t_a2 = calc_t.v(1.0 -
0.025)
L_Int = []
U_Int = []
for j in range(p):
L_Int.append(b[j, 0] -
sqrt(s2_b[j, j]) * t_a2)
U_Int.append(b[j, 0] +
sqrt(s2_b[j, j]) * t_a2)
f_out.write('\n')
f_out.write('95%
Confidence Intervals:\n')
for j in range(p):
f_out.write('b[{0}]<{1}>:\n
[{2:.5f}, {3:.5f}]\n'.\
format(j, XLabel[j], L_Int[j], U_Int[j]))
f_out.close()