第5章基本的回帰分析

スライド こちら

はじめに

これまでに第2章でデータ可視化、第3章でデータ操作、これらのスキルを基に、本章からはデータモデリングの世界へと進みます。データモデリングの基本的な考え方は、以下の関係性を明らかにすることです。

目的変数 (target variable) \(y\): 予測または説明したい変数。従属変数 (dependent variable) や 応答変数 (response variable) とも呼ばれます。
説明変数 (explanatory variable) / 予測変数 (predictor variable) \(x\): 目的変数を説明または予測するために用いる変数。独立変数 (independent variable) や 共変量 (covariate) とも呼ばれます。

数学的に言えば、目的変数 \(y\) を説明変数 \(x\) の「関数として」モデル化します。ここでの「関数」とは、数学的な意味での関数を指します。説明変数と予測変数という用語はしばしば同じ意味で使われますが、モデリングの目的によって使い分けられることがあります。データモデリングは、大まかに以下の二つの目的のいずれかのために行われます。

説明のためのモデリング (Modeling for Explanation): 目的変数 \(y\) と説明変数 \(x\) の間の関係性を明確に記述し、定量化することを目的とします。関係の統計的有意性を評価したり、変数間の関連性の強さを測る指標を得たり、場合によっては変数間の因果関係の理解を試みたりします。
予測のためのモデリング (Modeling for Prediction): 説明変数 \(x\) に含まれる情報に基づいて、目的変数 \(y\) の値を予測することを目的とします。この場合、変数間の関係性を深く理解することよりも、\(x\) を使って \(y\) を精度良く予測することに重点が置かれます。

本書では主に説明のためのモデリングに焦点を当て、変数を説明変数と呼びます。予測のためのモデリングに興味がある方は、『Pythonによる統計的学習入門』(An Introduction to Statistical Learning with Python) のような専門書を参照すると良いでしょう。世の中には決定木モデルやニューラルネットワークなど多種多様なモデリング手法が存在しますが、本書では最も基本的で理解しやすいアプローチの一つである線形回帰 (linear regression) に焦点を当てます。

線形回帰では、目的変数 \(y\) は数値型 (numerical) である必要があり、説明変数 \(x\) は数値型またはカテゴリ型 (categorical) のいずれかをとることができます。そして最も重要な仮定は、\(y\) と \(x\) の関係が線形、つまり直線的であると仮定することです。「直線」が何を意味するかは、説明変数 \(x\) の種類によって少し異なります。

本章（第5章）では、説明変数が1つだけのモデルを扱います。 * 5.1節では、説明変数が数値型である単純線形回帰 (simple linear regression) を学びます。 * 5.2節では、説明変数がカテゴリ型である場合を扱います。

続く第6章では、この基本的な回帰の考えを拡張し、2つの説明変数 \(x_1, x_2\) を持つモデルを検討します。そして第7章以降では、回帰モデルの結果を分析するために、第8章から第10章で学ぶ統計的推論 (statistical inference) の手法を用いていきます。

それでは、1つの説明変数 \(x\) を持つ線形回帰モデル、すなわち基本的な回帰分析から始めましょう。この過程で、相関係数 (correlation coefficient)、「相関は必ずしも因果を意味しない」という重要な注意点、そして「最適な当てはめ直線 (line of best fit)」とは何かといった統計学の基本概念についても解説していきます。

必要なパッケージ

この章で使用するPythonライブラリをインポートします。seabornライブラリから直接データセットを読み込むため、以前のようなデータダウンロードのステップは不要になります。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from IPython.display import display # DataFrameなどをきれいに表示するため

5.1 1つの数値的説明変数

レストランでは、食事の総支払額とウェイター/ウェイトレスへのチップの額にはどのような関係があるのでしょうか？総支払額が高いほど、チップの額も高くなる傾向があるのでしょうか？もしそうなら、その関係はどの程度の強さなのでしょうか？

この節では、レストランでの食事の総支払額という1つの数値型説明変数を用いて、チップの額の違いを説明することを試みます。これらの問いに答えるために、単純線形回帰を用いて、チップの額と食事の総支払額の関係をモデル化します。

目的変数 \(y\): チップの額 (tip)（数値型）
説明変数 \(x\): 食事の総支払額 (total_bill)（数値型）

使用するデータは、seabornライブラリに含まれるtipsデータセットです。このデータセットには、レストランのさまざまなテーブルでのチップの記録が含まれています。

探索的データ分析 (Exploratory Data Analysis, EDA)

まず、tipsデータセットを読み込み、基本的な情報を確認します。

# seabornからtipsデータセットを読み込む
tips = sns.load_dataset('tips')

EDAでは主に以下の3つのステップを行います。 1. データフレームの基本的な情報を確認し、実際のデータ値を見る。 2. 記述統計量（平均値、中央値、標準偏差など）を計算する。 3. データを可視化する。

1. 生のデータ値と構造の確認

# データの最初の数行を表示
print("レストランのチップデータの最初の5行:")
display(tips.head())

# データの構造（型や欠損値など）を確認
print("\nデータの構造:")
tips.info()

tipsデータセットには、total_bill (総支払額), tip (チップ額), sex (支払者の性別), smoker (喫煙席か否か), day (曜日), time (時間帯), size (グループの人数) といった列が含まれています。このセクションでは主に tip と total_bill に注目します。幸い、このデータセットには欠損値はないようです。

2. 記述統計量の計算

# tipとtotal_billの記述統計量の表示
print("\n'tip'と'total_bill'の記述統計量:")
display(tips[['tip', 'total_bill']].describe())

# チップ額と総支払額の相関係数を計算
correlation_tips = tips['tip'].corr(tips['total_bill'])
print(f"\nチップ額と総支払額の相関係数: {correlation_tips:.3f}")

相関係数は約0.676であり、チップ額と総支払額の間には比較的強い正の線形関係があることを示唆しています。つまり、総支払額が高いほど、チップの額も高い傾向があると言えそうです。

3. データの可視化

目的変数も説明変数も数値型なので、散布図を作成するのが適切です。

# 散布図の作成
plt.figure(figsize=(8, 5))
sns.scatterplot(x='total_bill', y='tip', data=tips, alpha=0.6)
plt.xlabel('総支払額 (total_bill)')
plt.ylabel('チップ額 (tip)')
plt.title('チップ額と総支払額の関係')
plt.grid(True, alpha=0.3)
plt.show()

この散布図から、点が右上がりの傾向をはっきりと示しており、これは先ほど計算した相関係数0.676（比較的強い正の相関）と整合的です。

回帰直線を重ねて全体の傾向を捉えやすくしてみましょう。

# 散布図と回帰直線の描画
plt.figure(figsize=(8, 5))
sns.regplot(x='total_bill', y='tip', data=tips,
            scatter_kws={'alpha':0.4, 's':50}, # 点の透明度とサイズ
            line_kws={'color':'blue'},        # 回帰直線の色
            ci=None)                          # 信頼区間は今回は表示しない
plt.xlabel('総支払額 (total_bill)')
plt.ylabel('チップ額 (tip)')
plt.title('チップ額と総支払額の関係 (回帰直線付き)')
plt.grid(True, alpha=0.3)
plt.show()

回帰直線は、データ点全体の右上がりの傾向をよく捉えています。

単純線形回帰

統計モデリング、特に回帰分析では、回帰直線の方程式は \(\widehat{y} = b_0 + b_1 x\) と書きます。 * \(b_0\) は切片 (intercept) で、\(x=0\) のときの \(\widehat{y}\) の値です。（ご質問のモデル \(y = a + bx\) における \(a\) に相当） * \(b_1\) は説明変数 \(x\) に対する傾き (slope) で、\(x\) が1単位増加したときの \(\widehat{y}\) の増加量です。（ご質問のモデル \(y = a + bx\) における \(b\) に相当） * \(\widehat{y}\) は、回帰モデルによって予測された予測値 (predicted value) または適合値 (fitted value) です。

Pythonのstatsmodelsライブラリを使って、線形回帰モデルを当てはめてみましょう。

# statsmodels を用いた線形回帰モデルの構築
# 'tip ~ total_bill' は「tip を total_bill で説明する」という意味の formula
model_tips_vs_bill = smf.ols('tip ~ total_bill', data=tips).fit()

# 回帰結果の要約を表示
print(model_tips_vs_bill.summary())

model_tips_vs_bill.summary() の出力結果の中の coef という列（係数の推定値）に注目します。 * Intercept (切片 \(b_0\)) の推定値は約 0.9203 です。 * total_bill (傾き \(b_1\)) の推定値は約 0.1050 です。

したがって、当てはめられた回帰直線の方程式は次のようになります。 \(\(\widehat{\text{tip}} = 0.9203 + 0.1050 \cdot \text{total\_bill}\)\)

係数の解釈:

切片 (\(b_0 \approx 0.9203\)): 「総支払額 total_bill」が0ドルの場合に予測されるチップ額が約0.92ドルであることを意味します。しかし、総支払額が0ドルというのは現実的な状況ではないため、この切片の値を実用的に解釈することにはあまり意味がありません。多くの場合、切片はモデルがデータに適合するために数学的に必要な値となります。
傾き (\(b_1 \approx 0.1050\)): 「総支払額 total_bill」が1ドル増加するごとに、チップ額 tip は平均して約0.105ドル（約10.5セント）増加する関連があると解釈されます。
- 符号が正であるため、これら2つの変数間には正の関係があることが示唆されます。これは相関係数 (0.676) が正であったことと整合します。
- 繰り返しになりますが、これはあくまで「関連」であり、この分析だけから「総支払額の増加がチップ額の増加を引き起こす」という因果関係を断定することはできません（5.3節で詳述）。
- 「平均して」という言葉が重要です。総支払額が1ドル増えたからといって、必ずチップが正確に0.105ドル増えるわけではなく、多くのデータで見た場合の平均的な傾向を示しています。

回帰結果の主要な情報をまとめた表（回帰表）を作成してみましょう。

# 回帰表の作成 (簡略版)
regression_table_tips_vs_bill = pd.DataFrame({
    'term': ['Intercept', 'total_bill'],
    'estimate': model_tips_vs_bill.params.values,
    'std_error': model_tips_vs_bill.bse.values,
    't_value': model_tips_vs_bill.tvalues.values,
    'p_value': model_tips_vs_bill.pvalues.values,
    'conf_int_lower': model_tips_vs_bill.conf_int().iloc[:, 0].values,
    'conf_int_upper': model_tips_vs_bill.conf_int().iloc[:, 1].values
})
print("\n回帰表 (tip ~ total_bill):")
display(regression_table_tips_vs_bill)

観測値、予測値、残差

回帰分析では、以下の3つの値が重要です。 1. 観測値 (\(y_i\)): 実際にデータとして観測された目的変数の値。 2. 予測値 (\(\widehat{y}_i\)): 当てはめられた回帰モデルから、ある説明変数の値 \(x_i\) に基づいて予測される目的変数の値。 3. 残差 (\(e_i = y_i - \widehat{y}_i\)): 観測値と予測値の差。モデルの「誤差」を表します。

例えば、最初のデータ点を見てみましょう。総支払額が \(x = 16.99\) ドル、実際のチップ額は \(y = 1.01\) ドルです。この場合の予測チップ額 \(\widehat{y}\) は、 \(\(\widehat{\text{tip}} = 0.9203 + 0.1050 \cdot 16.99 \approx 0.9203 + 1.7840 \approx 2.7043\)\)残差は、\(\(e = y - \widehat{y} = 1.01 - 2.7043 \approx -1.6943\)\) となります。この場合、モデルは実際のチップ額よりも高く予測したため、残差は負の値になっています。

これらの値を全てのデータ点について計算し、データフレームに追加してみましょう。

# 予測値と残差をデータフレームに追加
tips['tip_hat_vs_bill'] = model_tips_vs_bill.predict(tips)
tips['residual_vs_bill'] = model_tips_vs_bill.resid

# 最初の数行を表示して確認
print("\n予測値と残差を含むデータフレームの最初の数行:")
display(tips[['total_bill', 'tip', 'tip_hat_vs_bill', 'residual_vs_bill']].head())

残差プロット 残差の性質を調べるために、残差プロットを作成します。

# 残差プロットの作成
plt.figure(figsize=(8, 5))
sns.scatterplot(x='tip_hat_vs_bill', y='residual_vs_bill', data=tips, alpha=0.6)
plt.axhline(y=0, color='red', linestyle='--') # y=0 の線
plt.xlabel('予測チップ額 (ŷ)')
plt.ylabel('残差 (y - ŷ)')
plt.title('残差プロット (tip ~ total_bill)')
plt.grid(True, alpha=0.3)
plt.show()

理想的には、残差は \(y=0\) の線の周りにランダムに散らばり、明確なパターンが見られないことが望ましいです。このプロットでは、予測値が大きくなるにつれて残差のばらつきも若干大きくなる（扇形のような形）傾向が見られるかもしれません。これは不等分散性の兆候である可能性があり、より高度なモデリングでは考慮すべき点ですが、本章では深入りしません。

残差のヒストグラム

# 残差のヒストグラム
plt.figure(figsize=(8, 5))
sns.histplot(tips['residual_vs_bill'], kde=True, bins=20)
plt.xlabel('残差 (y - ŷ)')
plt.ylabel('度数')
plt.title('残差のヒストグラム (tip ~ total_bill)')
plt.grid(True, alpha=0.3)
plt.show()

残差は0を中心におおよそ対称的に分布しているように見えますが、少し右に裾が長い（正の方向に偏っている）かもしれません。

最適な当てはめ直線 (Line of Best Fit) と最小二乗法

回帰直線は「最適な当てはめ直線」とも呼ばれます。では、何をもって「最適」と判断するのでしょうか？回帰分析では、観測された全てのデータ点に対して、残差の平方和 (Sum of Squared Residuals, SSRまたはSSE) を最小にするような直線を「最適」な直線として選びます。残差平方和は以下のように定義されます。 \(\(SSR = \sum_{i=1}^{n} (y_i - \widehat{y}_i)^2 = \sum_{i=1}^{n} e_i^2\)\) ここで、\(n\) はデータ点の数です。各データ点での誤差（残差）を二乗し、それらを全て合計したものです。二乗することで、正の誤差も負の誤差も正の値として扱われ、大きな誤差ほどSSRへの寄与が大きくなります。このSSRを最小にするという基準は、最小二乗法 (Method of Least Squares) と呼ばれます。

最小二乗法による \(b_0\) と \(b_1\) の導出

単純線形回帰モデル \(\widehat{y}_i = b_0 + b_1 x_i\) において、SSRは次のように書けます。 \(\(SSR(b_0, b_1) = \sum_{i=1}^{n} (y_i - (b_0 + b_1 x_i))^2\)\) この \(SSR(b_0, b_1)\) を最小にする \(b_0\) と \(b_1\) を見つけるために、大学で学ぶ微分積分学の知識を使います。具体的には、\(SSR\) を \(b_0\) と \(b_1\) それぞれで偏微分し、その値を0とおいた連立方程式（正規方程式と呼ばれます）を解きます。

\(b_0\) に関する偏微分: \(\(\frac{\partial SSR}{\partial b_0} = \sum_{i=1}^{n} 2(y_i - b_0 - b_1 x_i)(-1) = -2 \sum_{i=1}^{n} (y_i - b_0 - b_1 x_i)\)\) これを0とおくと、 \(\(\sum_{i=1}^{n} (y_i - b_0 - b_1 x_i) = 0\)\) \(\(\sum y_i - n b_0 - b_1 \sum x_i = 0\)\) \(\(n b_0 + b_1 \sum x_i = \sum y_i \quad \cdots (1)\)\)
\(b_1\) に関する偏微分: \(\(\frac{\partial SSR}{\partial b_1} = \sum_{i=1}^{n} 2(y_i - b_0 - b_1 x_i)(-x_i) = -2 \sum_{i=1}^{n} x_i(y_i - b_0 - b_1 x_i)\)\) これを0とおくと、 \(\(\sum_{i=1}^{n} (x_i y_i - b_0 x_i - b_1 x_i^2) = 0\)\) \(\(b_0 \sum x_i + b_1 \sum x_i^2 = \sum x_i y_i \quad \cdots (2)\)\)

これで、\(b_0\) と \(b_1\) に関する2つの正規方程式が得られました。これらを解きます。まず、(1)式を \(n\) で割ると、 \(\(b_0 + b_1 \frac{\sum x_i}{n} = \frac{\sum y_i}{n}\)\)ここで、\(\bar{x} = \frac{\sum x_i}{n}\)（\(x\)の平均値）、\(\bar{y} = \frac{\sum y_i}{n}\)（\(y\)の平均値）なので、\(\(b_0 + b_1 \bar{x} = \bar{y}\)\) \(\(b_0 = \bar{y} - b_1 \bar{x} \quad \cdots (3)\)\) この(3)式は、回帰直線 \((\widehat{y} = b_0 + b_1 x)\) が必ず点 \((\bar{x}, \bar{y})\) を通ることを示しています。

次に、(3)式を(2)式に代入します。 \(\((\bar{y} - b_1 \bar{x}) \sum x_i + b_1 \sum x_i^2 = \sum x_i y_i\)\) \(\(\bar{y} \sum x_i - b_1 \bar{x} \sum x_i + b_1 \sum x_i^2 = \sum x_i y_i\)\)\(b_1\) について整理すると、\(\(b_1 (\sum x_i^2 - \bar{x} \sum x_i) = \sum x_i y_i - \bar{y} \sum x_i\)\)ここで、\(\bar{x} \sum x_i = \frac{(\sum x_i)^2}{n}\) および \(\bar{y} \sum x_i = \frac{(\sum y_i)(\sum x_i)}{n}\) を使うと、\(\(b_1 \left(\sum x_i^2 - \frac{(\sum x_i)^2}{n}\right) = \sum x_i y_i - \frac{(\sum y_i)(\sum x_i)}{n}\)\)分子・分母に \(n\) を掛けると、\(\(b_1 = \frac{n \sum x_i y_i - (\sum x_i)(\sum y_i)}{n \sum x_i^2 - (\sum x_i)^2}\)\)この式は、偏差平方和や偏差積和の記号 \(S_{xx} = \sum (x_i - \bar{x})^2 = \sum x_i^2 - n\bar{x}^2\) や \(S_{xy} = \sum (x_i - \bar{x})(y_i - \bar{y}) = \sum x_i y_i - n\bar{x}\bar{y}\) を用いると、より簡潔に以下のように書けます。\(\(b_1 = \frac{S_{xy}}{S_{xx}} = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2}\)\)そして、\(b_1\) が求まれば、(3)式から \(b_0\) が求まります。\(\(b_0 = \bar{y} - b_1 \bar{x}\)\)

これらの計算は複雑に見えるかもしれませんが、コンピュータ（例えばPythonのstatsmodelsライブラリ）が瞬時に行ってくれます。重要なのは、回帰直線がデータへの当てはまりの指標である残差平方和を最小にするように数学的に決定されている、という原理を理解することです。

# このモデルの残差平方和
ssr_tips_vs_bill = model_tips_vs_bill.ssr
print(f"このモデルの残差平方和 (SSR): {ssr_tips_vs_bill:.3f}")

statsmodelsが出力した係数 \(b_0\) と \(b_1\) は、まさにこのSSRを最小にする値です。

5.2 1つのカテゴリ的説明変数

レストランで支払われるチップの額は、曜日によって異なるのでしょうか？例えば、週末は平日よりもチップが多くなる傾向があるのでしょうか？

この節では、曜日という1つのカテゴリ型説明変数を用いて、チップの額の違いを探ります。引き続き tips データセットを使用します。

目的変数 \(y\): チップの額 (tip)（数値型）
説明変数 \(x\): 曜日 (day)（カテゴリ型: 木, 金, 土, 日）

探索的データ分析

1. day 変数の確認と記述統計量

print("\n曜日のカテゴリとそれぞれの件数:")
display(tips['day'].value_counts())

print("\n曜日ごとのチップ額の記述統計量:")
tip_by_day_summary = tips.groupby('day')['tip'].agg(['mean', 'median', 'std', 'min', 'max', 'count']).sort_values('mean', ascending=False)
display(tip_by_day_summary)

データには木曜日(Thur)、金曜日(Fri)、土曜日(Sat)、日曜日(Sun)の4つの曜日が含まれています。平均チップ額を見ると、日曜と土曜が高く、金曜が最も低いようです。

2. データの可視化

曜日ごとにチップ額の分布を箱ひげ図で比較します。

# 曜日別のチップ額の箱ひげ図
plt.figure(figsize=(8, 6))
# tips['day'].cat.categories でカテゴリの順序を取得して表示順を指定
day_order = tips['day'].cat.categories 
sns.boxplot(x='day', y='tip', data=tips, order=day_order)
plt.xlabel('曜日 (day)')
plt.ylabel('チップ額 (tip)')
plt.title('曜日別のチップ額の分布')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

# 曜日別のチップ額のバイオリンプロット
plt.figure(figsize=(8, 6))
sns.violinplot(x='day', y='tip', data=tips, order=day_order)
plt.xlabel('曜日 (day)')
plt.ylabel('チップ額 (tip)')
plt.title('曜日別のチップ額の分布 - バイオリンプロット')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

箱ひげ図やバイオリンプロットから、曜日によってチップ額の分布に違いが見られそうです。特に土曜日と日曜日は中央値が高く、ばらつきも大きいように見えます。

線形回帰

カテゴリ型の説明変数 day を用いて、チップ額 tip を説明する線形回帰モデルを当てはめてみましょう。

# カテゴリ変数を説明変数とした線形回帰モデル
model_tips_vs_day = smf.ols('tip ~ C(day)', data=tips).fit()

# 回帰結果の要約を表示
print(model_tips_vs_day.summary())

回帰結果の coef（係数）の部分に注目します。statsmodels は、tips データセットの day 列（カテゴリ型）のカテゴリの中で、通常、最初のもの（このデータセットでは Thur (木曜日)）を基準カテゴリ (baseline category) として自動的に選択します。 * Intercept: 約2.771。これは基準カテゴリである木曜日 (Thur) の平均チップ額の推定値です。 * C(day)[T.Fri]: 約-0.030。これは金曜日 (Fri) の平均チップ額が、木曜日の平均チップ額よりも平均して約0.030ドル低いことを示します。したがって、金曜日の平均チップ額の推定値は \(2.771 - 0.030 \approx 2.741\) ドルです。 * C(day)[T.Sat]: 約0.212。土曜日 (Sat) は木曜日より平均約0.212ドル高い。推定平均チップ額は \(2.771 + 0.212 \approx 2.983\) ドル。 * C(day)[T.Sun]: 約0.478。日曜日 (Sun) は木曜日より平均約0.478ドル高い。推定平均チップ額は \(2.771 + 0.478 \approx 3.249\) ドル。

他のカテゴリの係数は、この基準カテゴリとの差（オフセット）を表します。

回帰表を整理してみましょう。

# 回帰表の作成 (簡略版)
base_category = tips['day'].cat.categories[0] # 'Thur'
terms_day = [f'Intercept ({base_category})']
for term_name in model_tips_vs_day.params.index[1:]:
    day_name = term_name.split('[T.')[1][:-1]
    terms_day.append(f"{day_name} (vs {base_category})")

regression_table_tips_vs_day = pd.DataFrame({
    'term': terms_day,
    'estimate': model_tips_vs_day.params.values,
    'std_error': model_tips_vs_day.bse.values,
    't_value': model_tips_vs_day.tvalues.values,
    'p_value': model_tips_vs_day.pvalues.values,
    'conf_int_lower': model_tips_vs_day.conf_int().iloc[:, 0].values,
    'conf_int_upper': model_tips_vs_day.conf_int().iloc[:, 1].values
})
print(f"\nカテゴリ変数を用いた回帰の回帰表 (基準: {base_category}):")
display(regression_table_tips_vs_day)

この表の estimate を使うと、各曜日の平均チップ額の推定値が、EDAで計算した tip_by_day_summary の mean 列と一致することを確認できます。

観測値、予測値、残差

このモデルにおける予測値 \(\widehat{y}_i\) は、その食事がなされた曜日の平均チップ額の推定値となります。

# 予測値と残差をデータフレームに追加
tips['tip_hat_vs_day'] = model_tips_vs_day.predict(tips)
tips['residual_vs_day'] = model_tips_vs_day.resid

# 曜日、観測値、予測値、残差を表示 (最初の数行)
print("\n観測値、予測値、残差を含むデータフレーム (tip ~ day):")
display(tips[['day', 'tip', 'tip_hat_vs_day', 'residual_vs_day']].head())

例えば、最初のデータ点が日曜日 (Sun) で、実際のチップ額が1.01ドルだったとします。日曜日の予測平均チップ額 (tip_hat_vs_day) は約3.249ドルです。したがって、この場合の残差は約 \(1.01 - 3.249 = -2.239\) ドルとなります。

残差の分布を曜日ごとに見てみましょう。

# 曜日別の残差の箱ひげ図
plt.figure(figsize=(8, 6))
sns.boxplot(x='day', y='residual_vs_day', data=tips, order=day_order)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('曜日 (day)')
plt.ylabel('残差 (tip - tip_hat_vs_day)')
plt.title('曜日別の残差の分布')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

このプロットは、各曜日の平均チップ額からの個々のチップ額のばらつきを示しています。

5.3 相関は必ずしも因果を意味しない (Correlation does not imply causation)

この章を通じて、回帰係数を解釈する際には慎重な言葉遣いを心がけてきました。例えば、「総支払額 total_bill が1ドル増加するごとに、チップ額 tip は平均して約0.105ドル増加する関連がある」と述べました。「関連がある」という言葉を使うことで、直接的な因果関係を主張しているわけではないことを明確にしています。

2つの変数 X と Y の間に統計的な相関や関連が見られたとしても、それが必ずしも「X が Y を引き起こす」という因果関係を意味するわけではありません。この原則は非常に重要です。

例として、ある医師が「靴を履いたまま寝た患者は、翌朝頭痛を訴える傾向がある」という関連を医療記録から発見したとします。この医師が「靴を履いて寝ることが頭痛を引き起こす！」と結論付けるのは早計です。実際には、交絡変数 (confounding variable) が存在する可能性があります。この場合、それは「前夜の飲酒量」かもしれません。 * 飲酒量が多い (Z) \(\rightarrow\) 靴を履いたまま寝てしまう可能性が高い (X) * 飲酒量が多い (Z) \(\rightarrow\) 二日酔いで頭痛がする可能性が高い (Y)

この関係を図で示すと以下のようになります。

      Z (飲酒量)
     / \
    /   \
   v     v
  X       Y
(靴就寝) (頭痛)

この状況では、X と Y の間には見かけ上の関連がありますが、真の原因は Z である可能性があります。X と Y の真の因果関係を評価するためには、Z の影響を考慮に入れる（例：Z で調整する）必要があります。

因果関係を確立するには、ランダム化比較試験 (RCT) のような慎重に設計された実験的研究や、観察研究であっても高度な統計的手法（因果推論の手法）を用いて交絡変数の影響を制御する努力が必要です。線形回帰モデルを立てただけでは、通常、因果関係について強い主張をすることはできません。

まとめ

この章では、データモデリングの基本的なアプローチである線形回帰分析の基礎を学びました。特に、説明変数が1つの場合に焦点を当て、seabornライブラリのtipsデータセットを用いて具体的な分析を行いました。

主な学習ポイントは以下の通りです。 1. モデリングの目的: 「説明のためのモデリング」と「予測のためのモデリング」の違いを理解しました。 2. 単純線形回帰: 1つの数値型説明変数を用いて数値型の目的変数をモデル化する方法（\(\widehat{y} = b_0 + b_1 x\)）を学び、係数 \(b_0\)（切片）と \(b_1\)（傾き）の解釈を行いました。 3. 最小二乗法による係数の導出: 回帰直線の係数 \(b_0\) と \(b_1\) が、残差平方和を最小にするように数学的に導かれる過程を理解しました。 4. カテゴリ型説明変数を持つ回帰: 1つのカテゴリ型説明変数を用いて数値型の目的変数をモデル化する方法を学びました。基準カテゴリと比較したオフセットとして係数が解釈されることを見ました。 5. 探索的データ分析(EDA): モデリングの前段階として、データの可視化や記述統計量の確認が重要であることを再認識しました。 6. 観測値、予測値、残差: これらの基本的な概念と、残差プロットや残差の分布の確認の重要性を学びました。 7. 相関と因果: 「相関は必ずしも因果を意味しない」という統計学における重要な注意点を、交絡変数の例と共に学びました。

Pythonのライブラリ（pandas, matplotlib, seaborn, statsmodels）を使い、レストランのチップデータという身近な例でこれらの分析を実行し、結果を解釈する方法を体験しました。

次の章では、複数の説明変数を持つ重回帰分析 (multiple linear regression) へと進み、より複雑な関係性を捉えることができるモデルについて学んでいきます。

Pythonコード一括実行用 (Google Colab向け)

以下に、この章で用いたPythonコードをまとめて記述します。

# 必要なライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
from IPython.display import display

print("==================================================")
print("5.1 1つの数値的説明変数 (tipsデータセット)")
print("==================================================\n")

# seabornからtipsデータセットを読み込む
tips = sns.load_dataset('tips')

# 1. 生のデータ値と構造の確認
print("レストランのチップデータの最初の5行:")
display(tips.head())
print("\nデータの構造:")
tips.info()

# 2. 記述統計量の計算
print("\n'tip'と'total_bill'の記述統計量:")
display(tips[['tip', 'total_bill']].describe())
correlation_tips = tips['tip'].corr(tips['total_bill'])
print(f"\nチップ額と総支払額の相関係数: {correlation_tips:.3f}")

# 3. データの可視化
print("\n--- 散布図 ---")
plt.figure(figsize=(8, 5))
sns.scatterplot(x='total_bill', y='tip', data=tips, alpha=0.6)
plt.xlabel('総支払額 (total_bill)')
plt.ylabel('チップ額 (tip)')
plt.title('チップ額と総支払額の関係')
plt.grid(True, alpha=0.3)
plt.show()

print("\n--- 散布図と回帰直線 ---")
plt.figure(figsize=(8, 5))
sns.regplot(x='total_bill', y='tip', data=tips,
            scatter_kws={'alpha':0.4, 's':50},
            line_kws={'color':'blue'},
            ci=None)
plt.xlabel('総支払額 (total_bill)')
plt.ylabel('チップ額 (tip)')
plt.title('チップ額と総支払額の関係 (回帰直線付き)')
plt.grid(True, alpha=0.3)
plt.show()

# 単純線形回帰
print("\n--- 単純線形回帰モデル (tip ~ total_bill) ---")
model_tips_vs_bill = smf.ols('tip ~ total_bill', data=tips).fit()
print(model_tips_vs_bill.summary())

regression_table_tips_vs_bill = pd.DataFrame({
    'term': ['Intercept', 'total_bill'],
    'estimate': model_tips_vs_bill.params.values,
    'std_error': model_tips_vs_bill.bse.values,
    't_value': model_tips_vs_bill.tvalues.values,
    'p_value': model_tips_vs_bill.pvalues.values,
    'conf_int_lower': model_tips_vs_bill.conf_int().iloc[:, 0].values,
    'conf_int_upper': model_tips_vs_bill.conf_int().iloc[:, 1].values
})
print("\n回帰表 (tip ~ total_bill):")
display(regression_table_tips_vs_bill)

# 観測値、予測値、残差
tips['tip_hat_vs_bill'] = model_tips_vs_bill.predict(tips)
tips['residual_vs_bill'] = model_tips_vs_bill.resid

print("\n予測値と残差を含むデータフレームの最初の数行:")
display(tips[['total_bill', 'tip', 'tip_hat_vs_bill', 'residual_vs_bill']].head())

print("\n--- 残差プロット (tip ~ total_bill) ---")
plt.figure(figsize=(8, 5))
sns.scatterplot(x='tip_hat_vs_bill', y='residual_vs_bill', data=tips, alpha=0.6)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('予測チップ額 (ŷ)')
plt.ylabel('残差 (y - ŷ)')
plt.title('残差プロット (tip ~ total_bill)')
plt.grid(True, alpha=0.3)
plt.show()

print("\n--- 残差のヒストグラム (tip ~ total_bill) ---")
plt.figure(figsize=(8, 5))
sns.histplot(tips['residual_vs_bill'], kde=True, bins=20)
plt.xlabel('残差 (y - ŷ)')
plt.ylabel('度数')
plt.title('残差のヒストグラム (tip ~ total_bill)')
plt.grid(True, alpha=0.3)
plt.show()

ssr_tips_vs_bill = model_tips_vs_bill.ssr
print(f"\nこのモデルの残差平方和 (SSR): {ssr_tips_vs_bill:.3f}")


print("\n\n==================================================")
print("5.2 1つのカテゴリ的説明変数 (tipsデータセット)")
print("==================================================\n")

# 1. `day` 変数の確認と記述統計量
print("\n曜日のカテゴリとそれぞれの件数:")
display(tips['day'].value_counts())
print("\n曜日ごとのチップ額の記述統計量:")
tip_by_day_summary = tips.groupby('day')['tip'].agg(['mean', 'median', 'std', 'min', 'max', 'count']).sort_values('mean', ascending=False)
display(tip_by_day_summary)

# 2. データの可視化
day_order = tips['day'].cat.categories
print("\n--- 曜日別のチップ額の箱ひげ図 ---")
plt.figure(figsize=(8, 6))
sns.boxplot(x='day', y='tip', data=tips, order=day_order)
plt.xlabel('曜日 (day)')
plt.ylabel('チップ額 (tip)')
plt.title('曜日別のチップ額の分布')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

print("\n--- 曜日別のチップ額のバイオリンプロット ---")
plt.figure(figsize=(8, 6))
sns.violinplot(x='day', y='tip', data=tips, order=day_order)
plt.xlabel('曜日 (day)')
plt.ylabel('チップ額 (tip)')
plt.title('曜日別のチップ額の分布 - バイオリンプロット')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

# 線形回帰
print("\n--- カテゴリ変数を説明変数とした線形回帰モデル (tip ~ C(day)) ---")
model_tips_vs_day = smf.ols('tip ~ C(day)', data=tips).fit()
print(model_tips_vs_day.summary())

base_category = tips['day'].cat.categories[0]

terms_day = [f'Intercept ({base_category})']
for term_name in model_tips_vs_day.params.index[1:]:
    day_name = term_name.split('[T.')[1][:-1]
    terms_day.append(f"{day_name} (vs {base_category})")

regression_table_tips_vs_day = pd.DataFrame({
    'term': terms_day,
    'estimate': model_tips_vs_day.params.values,
    'std_error': model_tips_vs_day.bse.values,
    't_value': model_tips_vs_day.tvalues.values,
    'p_value': model_tips_vs_day.pvalues.values,
    'conf_int_lower': model_tips_vs_day.conf_int().iloc[:, 0].values,
    'conf_int_upper': model_tips_vs_day.conf_int().iloc[:, 1].values
})
print(f"\nカテゴリ変数を用いた回帰の回帰表 (基準: {base_category}):")
display(regression_table_tips_vs_day)

# 観測値、予測値、残差
tips['tip_hat_vs_day'] = model_tips_vs_day.predict(tips)
tips['residual_vs_day'] = model_tips_vs_day.resid

print("\n観測値、予測値、残差を含むデータフレーム (tip ~ day):")
display(tips[['day', 'tip', 'tip_hat_vs_day', 'residual_vs_day']].head())

print("\n--- 曜日別の残差の箱ひげ図 ---")
plt.figure(figsize=(8, 6))
sns.boxplot(x='day', y='residual_vs_day', data=tips, order=day_order)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('曜日 (day)')
plt.ylabel('残差 (tip - tip_hat_vs_day)')
plt.title('曜日別の残差の分布')
plt.grid(True, axis='y', alpha=0.3)
plt.show()

print("\n\n講義資料は以上です。")

第5章 基本的回帰分析

はじめに

必要なパッケージ

5.1 1つの数値的説明変数

探索的データ分析 (Exploratory Data Analysis, EDA)

単純線形回帰

観測値、予測値、残差

最適な当てはめ直線 (Line of Best Fit) と最小二乗法

5.2 1つのカテゴリ的説明変数

探索的データ分析

線形回帰

観測値、予測値、残差

5.3 相関は必ずしも因果を意味しない (Correlation does not imply causation)

まとめ

Pythonコード一括実行用 (Google Colab向け)

第5章基本的回帰分析