「入門はじめての統計解析」をPythonで写経 Vol.18 ～ ６章「はじめての回帰分析」①回帰直線の求め方

６章「はじめての回帰分析」

書籍の著者　石村貞夫 先生

この記事は、書籍「入門はじめての統計解析」６章「はじめての回帰分析」の Python写経活動 を取り扱います。

書籍の図・表・計算を淡々とPython化する写経シリーズです。
６章は 単回帰分析 に取り組みます。

この記事は 回帰直線をさまざま方法で計算 します！
ChatGPTの活用も継続してまいります！

では書籍を開いて統計解析の旅に出発です🚀

はじめに

このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

６章 はじめての回帰分析

この記事は６章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

６章で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd
import sympy

# 統計、回帰分析
import scipy.stats as stats
import pingouin as pg
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
import statsmodels.formula.api as smf
import lmdiag

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

イントロダクション

ChatGPTがストーリー～「線の物語」の語り部になります。

✅ 回帰分析の世界へようこそ

「この２つの変数、なんとなく関係がありそうだけど…具体的にどんな関係なんだろう？」

データを眺めていて、ふとそんな問いが浮かんだことはありませんか？
回帰分析は、２つの変数の“関係性のかたち”を、数式で表現してくれる手法です。
相関だけではつかみきれない「傾き」や「予測」のニュアンスも、回帰分析なら見えてきます。
今回から２回にわたり、回帰直線の求め方から、決定係数や分散分析表まで、基礎からしっかりと歩いていきます。

回帰分析の“はじめの一歩”を一緒に踏み出してみましょう。

➡️スタート：Section 6.1 回帰直線の求め方

テキストは２つの変数の間の関係を探る「単回帰分析」を取り扱います。
回帰直線は、２つの変数の関係を「直線」で表現するものです。
学校の数学で学んだ「切片」「傾き」が登場します。

ChatGPTが回帰分析の紹介文を作ってくれました。

✅ 回帰分析と回帰直線

２つの変数が関係していそうなとき、
「どんな関係なのかを、まっすぐな線で表してみよう」――
それが回帰分析のはじまりです。

商品価格と売上、気温とアイスの売れ行き、勉強時間とテストの点数……
“ある数値が、別の数値にどう影響しているのか” を知りたいとき、回帰分析はとても頼りになります。
さらに、その関係が式で表せれば、未来を予測することもできる。

実は回帰分析は、統計解析と機械学習が出会う場所でもあります。
データの傾向を理解し、予測に活かす。その最初の一歩として、回帰直線を一緒に体験してみましょう。

はじめてでも大丈夫。やさしく、のんびり進めていきます。

回帰直線のイメージ

テキスト p.248 表6.1.1 「宣伝広告費と売上高」をお借りして、データを読み込みます。

### データの登録 p.248 表6.1.1

# データフレーム
data1 = pd.DataFrame(
    {'会社名': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J'],
     '宣伝広告費': [107, 336, 233, 82, 61, 378, 129, 313, 142, 428],
     '売上高': [286, 851, 589, 389, 158, 1037, 463, 563, 372, 1020]},
     index=range(1, 11))
data1.index.name = 'No.'

# 結果の表示
print('data1.shape:', data1.shape)
data1

【実行結果】
10 社の宣伝広告費と売上高のデータです。

要約統計量を確認します。

### 要約統計量の表示
data1.describe().round(3)

【実行結果】
どの数値（countを除く）も売上高の方が大きいですね！

データを散布図で可視化して、２つの変数「宣伝広告費」「売上高」の関係を直感的に把握しましょう。
テキストの図 6.1.1「散布図」に相当します。
seaborn の scatterplot を利用します。

### 散布図の描画 p.248 図6.1.1
sns.scatterplot(data=data1, x='宣伝広告費', y='売上高', s=80)
plt.xlim(0, 520)
plt.ylim(0, 1250);

【実行結果】
横軸：宣伝広告費、縦軸：売上高のプロットです。
右上がりの傾向が見られます。
宣伝広告費が大きくなると売上高が大きくなる関係がありそうです。

相関係数を算出しましょう。
テキストの表 6.1.2「相関係数」に相当します。
pandas データフレームのメソッド corr() で相関行列を計算できます。

### 相関係数の算出 p.249 表6.1.2
data1.corr(numeric_only=True).round(3)

【実行結果】
宣伝広告費と売上高の間には強い正の相関があります！

では回帰直線 $${y=a + bx}$$ を描画しましょう。
テキストの図 6.1.2「回帰直線」に相当します。
ここで $${y}$$ は目的変数「売上高」、 $${x}$$ は説明変数「宣伝広告費」、$${a}$$ は切片、$${b}$$ は傾き（または回帰係数）です。

seaborn の lmplot() の機能で回帰直線を可視化します。
切片・傾きを計算しなくても回帰直線を描画できます！

### 回帰直線の描画 p.249 図6.1.2
sns.lmplot(data=data1, x='宣伝広告費', y='売上高', ci=None, height=4, aspect=1.5,
           scatter_kws={'s': 70}, line_kws={'color': 'tab:red', 'ls': '--'})
plt.text(x=180, y=750, s='$Y=a+bx$', fontsize=14)
plt.xlim(0, 520)
plt.ylim(0, 1250);

【実行結果】
赤い点線が回帰直線です。
青いデータ点の傍をほどよく通っていて、宣伝広告費と売上高の関係をうまく捉えているように感じます。

回帰直線を引くことができましたが…
せっかくなので切片・傾きを計算しましょう。

■ 回帰直線の傾きと切片の公式 p.251
テキストの公式をお借りします。
説明変数 $${x_i}$$、目的変数 $${y_i}$$、標本サイズ $${N}$$ とすると、切片 $${a}$$ と 傾き $${b}$$ は次のように計算できます。

$${\sum}$$ で埋め尽くされて圧迫感があります苦笑
よくよく式を見るとバリエーションは次の４つだけですのでご安心ください！

• 説明変数 $${x}$$ の合計 $${\sum_{i=1}^N x_i}$$

• 目的変数 $${y}$$ の合計 $${\sum_{i=1}^N y_i}$$

• 説明変数 $${x}$$ の二乗の合計 $${\sum_{i=1}^N x_i^2}$$

• 説明変数 $${x}$$ と目的変数 $${y}$$ の積合計 $${\sum_{i=1}^N x_i y_i}$$

【参考】公式・別バージョン
統計の書籍やWebサイトで見かける単回帰の公式の別バージョンです。

$${\bar{x}}$$：説明変数の標本平均
$${\bar{y}}$$：目的変数の標本平均
$${s_x}$$：説明変数の標本標準偏差
$${s_y}$$：目的変数の標本標準偏差
$${s_{xx}}$$：説明変数の標本分散
$${s_{xy}}$$：説明変数と目的変数の標本共分散
$${r_{xy}}$$：説明変数と目的変数の標本相関係数

回帰直線の求め方 例題 p.252

宣伝広告費と売上高のデータの切片と傾きを推定計算します。
さらに宣伝広告費が $${195}$$ のときの売上高を予測します。

説明変数が１つ、目的変数が１つの「単回帰」は、Python の多くのライブラリで実装できます。
主だったライブラリを使っていきましょう！

その前に、テキストの傾き・切片の公式を用いた関数を使ってみましょう。

🖲️自作関数
単回帰モデルの切片・傾きを算出する関数と予測関数を定義します。
４つの $${\sum}$$ の計算を「各種 Σ の算出」部分で行っています。

### 回帰直線の傾きと切片の公式

# 単回帰モデルの切片と傾きを算出する関数
def simple_linear_regerssion(X, Y):
    # 標本サイズの算出
    N = len(X)
    
    # 各種Σの算出
    sum_x = sum(X)
    sum_x2 = sum([x**2 for x in X])
    sum_y = sum(Y)
    sum_xy = sum([x * y for x, y in zip(X, Y)])
    
    # 傾きb, 切片aの算出
    b = (N * sum_xy - sum_x * sum_y) / (N * sum_x2 - sum_x**2)
    a = (sum_x2 * sum_y - sum_xy * sum_x) / (N * sum_x2 - sum_x**2)

    # 戻り値: 切片a, 傾きb
    return a, b

# 単回帰モデルの切片・傾きと説明変数xから目的変数yを予測する関数
def simple_linear_reg_prediction(x, intercept, slope):
    return intercept + slope * x

切片・傾きを推定します。

### 回帰直線の求め方 p.252
intercept, slope = simple_linear_regerssion(data1['宣伝広告費'], data1['売上高'])
intercept, slope

【実行結果】

回帰直線（回帰式）は次のようになります。

売上高の予測を行います。

# 売上高の予測

# 予測するxの値
x_new = 195

# 予測値の算出
pred = simple_linear_reg_prediction(x_new, intercept, slope)
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred:.2f} です')

【実行結果】
売上高の予測値は $${517}$$ です。

では Python のライブラリで回帰分析を実践しましょう。
お好みのライブラリが見つかるといいですね！

🖲️scipy.stats の linregress()
統計用のライブラリです。
切片・傾きを推定します。
引数は説明変数、目的変数です。

# scipy.stats利用
result1_stats = stats.linregress(data1['宣伝広告費'], data1['売上高'])
result1_stats[1], result1_stats[0]

【実行結果】
自作関数の結果と一致しました！

売上高の予測をします。
自作関数を使います。

# 売上高の予測
pred = simple_linear_reg_prediction(x_new, result1_stats[1], result1_stats[0])
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred:.2f} です')

【実行結果】
自作関数の結果と一致しました！

🖲️scikit-learn の LinearRegression()
機械学習のライブラリです。
切片・傾きを推定します。
引数は X=説明変数、y=目的変数です。

# scikit-learn利用
reg1 = LinearRegression()
reg1.fit(X=data1[['宣伝広告費']], y=data1['売上高'])
reg1.intercept_, reg1.coef_[0]

【実行結果】
自作関数の結果と一致しました！

売上高の予測をします。
学習済みモデル reg1 に対して predict を適用します。

# 売上高の予測
pred = reg1.predict(pd.DataFrame(dict(宣伝広告費=[x_new])))
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred[0]:.2f} です')

【実行結果】
自作関数の結果と一致しました！

🖲️pingouin の linear_regression()
統計用のライブラリです。
切片・傾きを推定します。
引数は X=説明変数、y=目的変数です。

# pingouin利用
result1_pg = pg.linear_regression(X=data1[['宣伝広告費']], y=data1['売上高'])
result1_pg

【実行結果】
自作関数の結果と一致しました！
さまざまな統計量を確認できます。

売上高の予測をします。
自作関数を使います。

# 売上高の予測
pred = simple_linear_reg_prediction(
    x_new, result1_pg.iloc[0, 1], result1_pg.iloc[1, 1])
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred:.2f} です')

【実行結果】
自作関数の結果と一致しました！

🖲️numpy の lstsq()
数値計算用ライブラリの最小二乗法関数です。
切片・傾きを推定します。
引数は X=説明変数、y=目的変数などです。
最初のコードでは、説明変数に定数項（すべての値が$${1}$$）を追加しています。

# numpy利用

# 説明変数Xに定数項を追加
X = np.vstack([np.ones(len(data1)), data1['宣伝広告費']]).T

# 最小二乗法の実行
result1_np = np.linalg.lstsq(X, data1['売上高'].values, rcond=-1)[0]
result1_np

【実行結果】
自作関数の結果と一致しました！

売上高の予測をします。
自作関数を使います。

# 売上高の予測
pred = simple_linear_reg_prediction(x_new, *result1_np)
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred:.2f} です')

【実行結果】
自作関数の結果と一致しました！

🖲️statsmodels の ols()
統計用のライブラリです。
R 言語風の formula API を使います。
切片・傾きを推定します。
引数は 回帰式を示す formula、データフレームです。

# statsmodels利用
result1_sm = smf.ols(formula='売上高 ~ 宣伝広告費', data=data1).fit()
result1_sm.summary()

【実行結果】
中段下の 「coef」列に 切片「Intercept」、傾き「宣伝広告費」が表示されています。
自作関数の結果と一致しました！
さまざまな統計量を確認できます。

売上高の予測をします。
学習済みモデル result1_sm に対して predict を適用します。

# 売上高を予測
pred = result1_sm.predict(dict(宣伝広告費=x_new))
print(f'宣伝広告費 {x_new} のときの売上高の予測値は {pred[0]:.2f} です')

【実行結果】
自作関数の結果と一致しました！

■ 残差の可視化
回帰モデルの残差は、目的変数の観測値（データの値）とモデルによる予測値の差です。
回帰モデルが満たすべき「仮定」をチェックするときに残差を確認します。
statsmodels の学習済みモデル を使って、lmdiag ライブラリで可視化することができます。

#　残差プロット 
# 参考サイト https://py4etrics.github.io/10_Residuals.html

plt.figure(figsize=(10, 6))
lmdiag.plot(result1_sm);

【実行結果】
概ね良好な感触です。

左の縦２つのチャートは残差がランダムになっているか（何か傾向がないか）を確認します。
右上のチャートは Vol.13 の正規性の検定で紹介した 正規 Q-Q プロットです。
残差が正規分布に従っているかを確認します。
右下のチャートは外れ値の存在を確認するチャートです。

■ 最小二乗法を計算する
回帰直線の切片・傾きを計算する際、残差の二乗和が最小になるように「最小二乗法」という計算方法を採っています。
最小二乗法の計算を Python で追いかけてみましょう。
「残差平方和 ▶️ 偏微分 ▶️ 連立方程式」の手順で進めます。

計算が大変なので、sympy ライブラリに計算を肩代わりしてもらいます。
こちらのWebサイトの情報を参考にいたしました。
ありがとうございます！

sympy の設定、残差平方和の算出を行います。

# 最小二乗法を１つ１つ計算する sympy利用
# 参考サイト: https://note.crohaco.net/2018/least-square-method/

# TeX表示の開始
sympy.init_printing()

# sympyの変数の定義 a:切片、b:傾き
a, b = sympy.symbols('a b')

# 残差平方和の算出
rss = sum([(y - (a + b*x))**2 for x, y in data1[['宣伝広告費', '売上高']].values])
rss.expand()

【実行結果】
テキスト p.251 の「残差の２乗和」を展開すると、この式になります。

残差平方和について、$${a}$$ と $${b}$$ の偏微分を求めます。
sympy.diff(関数, 変数) で変数の偏微分を計算できます。

# 残差平方和を切片a, 傾きbで偏微分
expr_a = sympy.diff(rss, a)
expr_b = sympy.diff(rss, b)
expr_a, expr_b

【実行結果】

偏微分$${=0}$$ とおいて連立方程式とし、連立方程式の解を求めます。
sympy.solve([方程式１, 方程式２]) で連立方程式を解きます。

# 偏微分を０として連立方程式を解く
solved = sympy.solve([expr_a, expr_b])
solved

【実行結果】
分数の形式で切片 $${a}$$、傾き $${b}$$ の解を得ました。

切片 $${a}$$、傾き $${b}$$ を小数の形式に変形します。

# 切片a, 傾きbの推定値の表示
[(sol[0], float(sol[1])) for sol in solved.items()]

【実行結果】
切片 $${a}$$、傾き $${b}$$ の値は自作関数の結果と一致しました。

TeX表示を解除します。

# TeX表示の解除
sympy.init_printing(pretty_print=False)

■ 分散共分散行列から傾き $${b}$$ を計算する
テキスト p.252 の右下で袴の人が語っています。

この公式は「公式・別バージョン」にも現れています！

同ページ上部には、ウサギの人が分散共分散行列の計算結果を示しています。
さあ分散共分散行列から、切片と傾きを計算しましょう。

分散共分散行列を計算します。
pandas データフレームに対して cov() を適用することでデータの分散共分散行列を作成できます。

### 分散共分散行列の算出 p.252
cov_matrix = data1[['宣伝広告費', '売上高']].cov(ddof=1)
cov_matrix.round(3)

【実行結果】
ウサギの人が語る分散共分散行列と一致しました！

続いて分散共分散行列から切片と傾きを計算します。
傾きは「共分散 ÷ 説明変数の分散」で計算できます。
切片は「目的変数の標本平均 － 傾きの推定値 × 説明変数の標本平均」で計算できます。

### 分散共分散行列で傾きを、x,yの平均値を用いて切片を推定 p.252

## 傾きb_hatの算出
# 分散共分散行列からxyの共分散とxの分散を取得
cov_xy, var_x = cov_matrix.iloc[0, 1], cov_matrix.iloc[0, 0]
# 傾きの算出 ※b_hat = xyの共分散 / xの分散
b_hat = cov_xy / var_x

## 切片の算出 ※a_hat = y_bar - b_hat * x_bar
# 標本平均x_bat, y_barの算出
x_bar, y_bar = data1[['宣伝広告費', '売上高']].mean()
# 切片の算出
a_hat = y_bar - b_hat * x_bar

## 結果の表示
a_hat, b_hat

【実行結果】
左：切片、右：傾きを計算できました。

理解度チェック 回帰直線の求め方 p.253

海水温と海洋生物の活動時間の関係を回帰分析で確かめます。

コードを淡々と書きます。
テキストのデータをお借りします。

### 回帰直線の求め方 海水の温度とクルマエビの夜間活動時間 p.253 表6.1.6

# データの登録
data2 = pd.DataFrame(
    {'No.': range(1, 11),
     '海水温度': [30.4, 27.2, 30.9, 22.5, 19.0, 16.4, 12.1, 12.7, 13.7, 23.6],
     '活動時間': [5.7, 6.7, 7.6, 7.7, 6.9, 4.6, 3.6, 6.4, 7.5, 6.4]})

# 結果の表示
print('data2.shape: ', data2.shape)
data2

【実行結果】
データの個数（標本サイズ）は $${10}$$ です。

散布図で２つの変数の関係を直感的に把握しましょう。

# 散布図で可視化
plt.figure(figsize=(7, 4))
sns.scatterplot(data=data2, x='海水温度', y='活動時間', s=70)
plt.xlim(10, 34)
plt.ylim(0, 10);

【実行結果】
目を細めると「水平に近い回帰直線」が見えてきませんか…

相関係数を確認します。

# 相関係数の確認
data2[['海水温度', '活動時間']].corr().round(3)

【実行結果】
相関係数 $${0.384}$$ は弱い正の相関です。

それでは回帰直線を求めましょう。
自作関数で切片と傾きの推定値を求めて、回帰式を表示します。

### 単回帰分析の実行 p.253

# 関数利用
a, b = simple_linear_regerssion(data2['海水温度'], data2['活動時間'])
print(f'回帰直線 Y = {a:.3f} {b:+.3f}x')

【実行結果】
傾きは小さな値です。

データの散布図に回帰直線を重ねて描画しましょう。

# 散布図と回帰直線を可視化

# 単回帰を用いた予測を実行する関数の定義
f = lambda x: a + b * x

# 描画領域の設定
plt.figure(figsize=(7, 4))
# 観測値の散布図の描画
sns.scatterplot(data=data2, x='海水温度', y='活動時間', s=70)
# 回帰直線の描画
plt.plot([10, 34], [f(10), f(34)], color='tab:red', ls='--')
# 修飾
plt.xlim(10, 34)
plt.ylim(0, 10);

【実行結果】
海水温度が高くなると活動時間が長くなるような関係が見えました。

ここからは Python のさまざまなライブラリで切片・傾きの推定値を軽快に求めていきましょう♪

🖲️scipy.stats

# scipy.stats利用
result2_stats = stats.linregress(data2['海水温度'], data2['活動時間'])
print('切片:', result2_stats[1], ', 傾き:', result2_stats[0])

【実行結果】

🖲️scikit-learn

# scikit-learn利用
reg2 = LinearRegression()
reg2.fit(data2[['海水温度']], data2['活動時間'])
print('切片:', reg2.intercept_, ', 傾き:', reg2.coef_[0])

【実行結果】

🖲️pingouin

# pingouin利用
result2_pg = pg.linear_regression(data2[['海水温度']], data2['活動時間'])
result2_pg

【実行結果】

🖲️numpy

# numpy利用

# 説明変数Xに定数項を追加
X = np.vstack([np.ones(len(data2)), data2['海水温度']]).T
# 最小二乗法の実行
a, b = np.linalg.lstsq(X, data2['活動時間'].values, rcond=-1)[0]
print('切片:', a, ', 傾き:', b)

【実行結果】

🖲️statsmodels

# statsmodels利用
result2_sm = smf.ols(formula='活動時間 ~ 海水温度', data=data2).fit()
a, b = result2_sm.params
print('切片:', a, ', 傾き:', b)

【実行結果】

残差プロットも描画しておきましょう。

#　残差プロット 
plt.figure(figsize=(10, 6))
lmdiag.plot(result2_sm);

【実行結果】
特段、問題となるようなことは見られません。

記事の最後をChatGPTに締めくくってもらいましょう！
「線の物語」の一部を綴っています。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。