単回帰分析
単回帰分析を行うプログラムを作成した。計算式は、＜岡本安晴「心理学データ分析と測定」勁草書房＞の第１章におけるものを用いた。
上掲書の表1.1.1のデータに対して、単回帰モデル
を当てはめる。
|
表1.1.１ 英語と国語の点数。（英語、国語）の形式で表す。 |
|
(60, 71) (61, 53) (48, 61) (55, 57) (66, 58) (62, 65) (50, 48) (51, 56) (57, 59) (66, 70) (71, 67) (71, 73) (67, 54) (76, 72) (67, 66) (60, 66) (48, 54) (75, 70) (61, 60) (60, 62) (58, 64) (75, 76) (68, 65) (66, 69) (67, 62) (65, 60) (74, 70) (72, 74) (38, 39) (65, 71) (50, 49) (51, 51) (53, 53) (64, 54) (68, 62) (51, 51) (61, 64) (74, 66) (58, 61) (28, 40) (78, 72) (64, 68) (62, 63) (52, 53) (70, 70) (78, 79) (75, 74) (59, 65) (73, 76) (63, 62) |
上のデータに単回帰モデル（１）を最小二乗法の基準で適用して分析するメイン プログラムのリストは、本ウェブサイトの後半に示した。プログラムのファイルはまとめて圧縮ファイルSmplRgrsnPrg.zipとした。表1.1.1のデータは、図１の形式で用意する。
・
・
・
図１
第１文字目がスラッシュ / で始まる行の次の行に、従属変数名（図１ではEng）と独立変数名（図１ではLng）をタブ文字あるいは半角空白文字で区切って並べて書く。第１文字目がスラッシュ / で始まる行の前までは、ファイルは読み飛ばされるので、コメントなどを書くことができる。
入力データファイルの準備後、プログラムを実行する。入力データファイル名の設定が求められる（図２）。
図２
プログラムと同じフォルダ内にデータファイルが置かれている場合は、ファイル名のみの指定でよい（図３）。
図３
入力データファイル名を設定して、Enterキーを押すと、次に出力ファイル名の設定が求められる（図４）。
図４
出力ファイル名を設定して、Enterキーを押すと、計算が始まり、図５のグラフが表示される。
図５
図５のウィンドウの右上のＸをクリックして閉じると、プログラムの実行終了となる。
プログラムの実行終了後、出力ファイルを開くと、表1.1.1のデータの場合、以下のような内容になっている。
Input data file =
Data111.txt
Input data...
60.0 71.0
61.0 53.0
48.0 61.0
55.0 57.0
66.0 58.0
・
・
・
70.0 70.0
78.0 79.0
75.0 74.0
59.0 65.0
73.0 76.0
63.0 62.0
b0 =
1.2636356886523998
b1 =
0.9756218289815616
Eng = 1.263636 +
0.9756218 * Lng
Coefficient of
determination: R2 = 0.7318389
Correlation
coefficient:
r = 0.8554758
t = 11.44539 p-value = 0.0000000 df = 48
95% Confidence
Interval = [0.8042323, 1.147011]
１行目に、読み込んだデータファイル名が出力されている。
次に、読み込んだデータが出力されている。これにより、データが正しく読み込まれたかどうかのチェックを行うことができる。
データの出力後、回帰モデル（１）の係数の値、b0（回帰直線のＹ切片）とb1（傾き）、が出力されている。これらの値を回帰モデル（１）に代入したときの回帰式も
Eng = 1.263636 + 0.9756218 * Lng
と出力されている。
続いて、決定係数と相関係数が
Coefficient of determination: R2 =
0.7318389
Correlation coefficient: r = 0.8554758
と出力されている。
次に、係数b1に関するｔの値と、帰無仮説「b1=0」に対するｐ値（両側検定）が
t = 11.44539 p-value = 0.0000000 df = 48
と出力されている。
最後に、９５％信頼区間が
95% Confidence Interval = [0.8042323, 1.147011]
と出力されている。
以下に、プログラムのメインソースファイルを示す。メインソースファイルを含めたプログラムファイル、およびデータファイルData111.txtは圧縮ファイルとしてまとめられている。
from math import *
from splitstring
import *
from calc_p_from_t
import *
from calc_t_from_p
import *
import matplotlib
import
matplotlib.pyplot as plt
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')
splitStr.set_string(data_f[pos])
YLabel =
splitStr.seek_nonempty(0)
XLabel =
splitStr.seek_nonempty(1)
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(float(splitStr.seek_nonempty(1)))
n += 1
pos += 1
f_out.write("
")
f_out.write("Input
data...\n")
sumX = 0.0
sumY = 0.0
for i in range(n):
f_out.write("{0:10}{1:10}\n".format(Y[i], X[i]))
sumX += X[i]
sumY += Y[i]
meanX = sumX / n
meanY = sumY / n
SSX = 0.0
SSY = 0.0
SXY = 0.0
for i in range(n):
SSX += (X[i] - meanX) ** 2.0
SSY += (Y[i] - meanY) ** 2.0
SXY += (X[i] - meanX) *
(Y[i] - meanY)
b1 = SXY / SSX
b0 = meanY - b1 *
meanX
f_out.write("\nb0
= {}\n".format(b0))
f_out.write("b1 =
{}\n".format(b1))
f_out.write("{0}
= {1:.7} + {2:.7} * {3}\n".format(YLabel, b0, b1, XLabel))
r = SXY / (sqrt(SSX) * sqrt(SSY))
R2 = r ** 2.0
f_out.write("\nCoefficient
of determination: R2 = {0:.7}\n".format(R2))
f_out.write("Correlation
coefficient:
r = {0:.7}\n".format(r))
SSE = 0.0
for i in range(n):
SSE += (Y[i] - (b0 + b1 *
X[i])) ** 2.0
MSE = SSE / (n - 2.0)
s2b1 = MSE / SSX
sb1 = sqrt(s2b1)
t = b1 / sb1
calc_p = C_cum_t(n -
2)
p = calc_p.v(t)
if p < 0.5:
p = 2 * p
else:
p = 2 * (1.0 - p)
f_out.write("\nt
= {0:.7} p-value =
{1:.7f} df = {2}\n".\
format(t, p, n - 2))
calc_t = CCalc_t(n -
2)
t_a2 = calc_t.v(1.0 -
0.025)
L_Int = b1 - sb1 *
t_a2
U_Int = b1 + sb1 *
t_a2
f_out.write("\n95%
Confidence Interval = [{0:.7}, {1:.7}]\n".format(L_Int, U_Int))
f_out.close()
min_x = X[0]
max_x = X[0]
for i in range(n):
if min_x > X[i]:
min_x = X[i]
if max_x < X[i]:
max_x = X[i]
L_y = b0 + b1 * min_x
R_y = b0 + b1 * max_x
plt.xlabel(XLabel)
plt.ylabel(YLabel)
plt.title("{0} =
{1:.7} + {2:.7} * {3}".format(YLabel, b0, b1, XLabel) +\
"R2 = {0:.5} r =
{1:.5}".format(R2, r))
plt.plot(X, Y, 'go')
plt.plot([min_x,
max_x], [L_y, R_y])
plt.show()