第１章 ２つの変数の間の関係を見る
単回帰分析
2つの変量の大まかな関係が1次式で表されることがあり、単回帰モデルと呼ばれている。図１のデータ（ファイル名：data1_1.xlsx）を考える。
図１
図１のデータは、＜岡本安晴「心理学データ分析と測定」勁草書房＞の表1.1.1の一部を抜き出したものである。図１のデータは、次のスクリプトで読み込むことができる。
xlsx =
pd.ExcelFile('data1_1.xlsx') # 入力ファイル名の設定
data = pd.read_excel(xlsx) # 読み込み
読み込まれたデータはDataFrame型としてdataで表されている。このdataからExcelファイルの１行目に設定されている変数名と、それぞれの変数名の列に設定されているデータ値を読み込み、散布図を描くスクリプトは以下のようになる。
CaseLabel = data.columns[0]
# 第０列目の変数名
CaseID = data[CaseLabel].values
# CaseLabelの表す名前の列のデータを取り出してCaseIDで表す
print(CaseID)
YLabel = data.columns[1] # 第１列目の変数名
Y = data[YLabel].values
# YLabelの表す名前の列のデータを取り出してYで表す
print(Y)
XLabel = data.columns[2]
# 第２列目の変数名
X = data[XLabel].values
# XL+abelの表す名前の列のデータを取り出してXで表す
print(X)
#
#
(X, Y)の散布図
#
plt.plot(X, Y, 'bo')
plt.xlabel(XLabel)
plt.ylabel(YLabel)
plt.show()
上のスクリプトの実行により図２の散布図が描かれる。
図２
変量Xの値が増加すると変量Yの値も増加する傾向にあることがわかる。この傾向を次式の単回帰モデルで表す。
式番号（1.3.1）は＜岡本安晴「心理学データ分析と測定」勁草書房＞におけるものである。本ウェブサイトにおける式番号などは、本拙著におけるものを用いる。
回帰モデルの係数は次式
で与えられるが、この計算は次のスクリプトで行っている。
X_mean = X.mean()
Y_mean = Y.mean()
num = ((X - X_mean) * (Y -
Y_mean)).sum()
den = ((X - X_mean) ** 2).sum()
b1 = num / den
b0 = Y_mean - b1 * X_mean
求められた回帰係数による回帰直線は、次のスクリプトで描いている。
plt.plot(X, Y, 'bo')
plt.plot([X.min(), X.max()], [b0 + b1
* X.min(), b0 + b1 * X.max()], 'g')
plt.title('b0 = {0:.3}, b1 = {1:.3}'.format(b0, b1), fontsize =
16)
plt.xlabel(XLabel, fontsize = 14)
plt.ylabel(YLabel, fontsize = 14)
plt.show()
図２のウィンドウの右上のXアイコンをクリックして閉じ、上の回帰直線を描くスクリプトを実行すると図３のグラフが描かれる。
図３
描かれた図は、「保存」ボタンのクリックでファイルに保存できる。
決定係数は次式で与えられる。
これは、次のスクリプトで計算できる。
SSReg = (((b0 + b1 * X) - Y_mean) **
2).sum()
SSTol = ((Y - Y_mean) ** 2).sum()
R2 = SSReg / SSTol
帰無仮説「
」に対するｔ値と、その両側検定におけるｐ値は、次のスクリプトで計算できる。
n = len(Y)
SSE = ((Y - (b0 + b1 * X)) ** 2).sum()
MSE = SSE / (n - 2)
SSX = ((X - X_mean) ** 2).sum()
s2_b1 = MSE / SSX
s_b1 = s2_b1 ** 0.5
t = b1 / s_b1
pvalue = (1 - ss.t.cdf(t, n - 2)) *
2
信頼区間（式（1.3.8））を求めるスクリプトは、以下のようである。
t_975 = ss.t.ppf(0.975, n - 2)
print('\n95% CI = [{0:<.5},
{1:<.5}]'.format(b1 - t_975 * s_b1, b1 + t_975 * s_b1))
上のスクリプトを1つにまとめたものをリスト１に示す。ファイルは、圧縮ファイルfiles1.zipにまとめたので、ダウンロードすれば自由に利用できる。
リスト１ 単回帰分析のPythonスクリプト例
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
xlsx =
pd.ExcelFile('data1_1.xlsx') # 入力ファイル名の設定
data = pd.read_excel(xlsx)
# 読み込み
CaseLabel = data.columns[0]
# 第０列目の変数名
CaseID = data[CaseLabel].values
# CaseLabelの表す名前の列のデータを取り出してCaseIDで表す
print(CaseID)
YLabel = data.columns[1]
# 第１列目の変数名
Y = data[YLabel].values
# YLabelの表す名前の列のデータを取り出してYで表す
print(Y)
XLabel = data.columns[2]
# 第２列目の変数名
X = data[XLabel].values
# XL+abelの表す名前の列のデータを取り出してXで表す
print(X)
#
#
(X, Y)の散布図
#
plt.plot(X, Y, 'bo')
plt.xlabel(XLabel)
plt.ylabel(YLabel)
plt.show()
#
# b0とb1の計算
#
X_mean = X.mean()
Y_mean = Y.mean()
num = ((X - X_mean) * (Y -
Y_mean)).sum()
den = ((X - X_mean) ** 2).sum()
b1 = num / den
b0 = Y_mean - b1 * X_mean
print('\nb0 = ', b0)
print('b1 = ', b1)
#
# 回帰直線と散布図の描画
#
plt.plot(X, Y, 'bo')
plt.plot([X.min(), X.max()], [b0 + b1
* X.min(), b0 + b1 * X.max()], 'g')
plt.title('b0 = {0:.3}, b1 = {1:.3}'.format(b0, b1), fontsize =
16)
plt.xlabel(XLabel, fontsize = 14)
plt.ylabel(YLabel, fontsize = 14)
plt.show()
#
# 決定係数R^2の計算
#
SSReg = (((b0 + b1 * X) - Y_mean) **
2).sum()
SSTol = ((Y - Y_mean) ** 2).sum()
R2 = SSReg / SSTol
print('\nR2 = ', R2)
#
# t値の計算
#
n = len(Y)
SSE = ((Y - (b0 + b1 * X)) ** 2).sum()
MSE = SSE / (n - 2)
SSX = ((X - X_mean) ** 2).sum()
s2_b1 = MSE / SSX
s_b1 = s2_b1 ** 0.5
t = b1 / s_b1
pvalue = (1 - ss.t.cdf(t, n - 2)) *
2
# 両側検定
print('t = ', t, ' p = ', pvalue)
#
# 95%信頼区間
#
t_975 = ss.t.ppf(0.975, n - 2)
print('\n95% CI = [{0:<.5},
{1:<.5}]'.format(b1 - t_975 * s_b1, b1 + t_975 * s_b1))