5. Multiple Regression

第6章重回帰分析 (Multiple Regression Analysis)

学習目標

この章を学ぶことで、以下の能力を習得します。

複数の説明変数を用いて目的変数を予測する重回帰モデルを定式化し、解釈する。
数値変数とカテゴリ変数が混在するデータに対して、交互作用モデルと平行傾斜モデルを構築し、その違いを理解する。
回帰平面の概念を理解し、2つの数値説明変数を持つモデルを可視化・解釈する。
モデル選択の基本的な考え方を理解し、可視化や決定係数 \(R^2\) を用いてモデルの適合度を評価する。
シンプソンのパラドックスのような、多変量データに潜む統計的現象を特定し、説明する。

スライド こちら データセット こちら Kaggleリンク こちら

導入

第5章では、単一の説明変数を用いて結果をモデル化する単純線形回帰を扱いました。しかし、現実世界の現象の多くは、単一の要因だけで説明できるほど単純ではありません。例えば、レストランのチップの額は、会計総額だけでなく、サービスを受けた客の性別、曜日、時間帯など、複数の要因に影響されると考えるのが自然です。

本章では、この考え方を拡張し、重回帰分析 (Multiple Regression Analysis) の世界へと踏み出します。重回帰は、複数の説明変数 \(x\_1, x\_2, \\dots, x\_p\) を同時に考慮に入れて目的変数 \(y\) をモデル化する強力な手法です。

複数の変数を同時に分析することで、各変数が他に与える影響を調整した上での「純粋な」関連性を評価できます。これは、しばしば交絡（confounding）によって引き起こされる見せかけの相関を排し、より本質的な関係性を探る上で不可欠です。本章を通じて、多面的なデータから知見を引き出すための分析スキルを構築していきましょう。

必要なライブラリ

本章の分析では、以下のPythonライブラリを使用します。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression
# 決定係数や平均二乗誤差の計算に使用
from sklearn.metrics import r2_score, mean_squared_error
# 3Dプロット用
from mpl_toolkits.mplot3d import Axes3D

6.1 1つの数値変数と1つのカテゴリ変数を持つモデル

ここでは、seabornライブラリに付属するtipsデータセットを用います。このデータは、レストランでのチップの支払いに関する情報を含んでいます。

今回は、以下の説明変数の組み合わせでチップ額をモデル化します。

目的変数 \(y\) (数値): チップ額 (tip)
説明変数 (2つ):
1. \(x\_1\) (数値): 会計総額 (total_bill)
2. \(x\_2\) (カテゴリ): 支払い客の性別 (sex)

会計総額はチップ額にどう影響するでしょうか？また、その影響は性別によって異なるのでしょうか？これらの問いに答えるため、重回帰モデルを適用します。

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

まず、データを読み込み、分析の準備を整えます。

# seabornからtipsデータセットをロード
tips = sns.load_dataset("tips")

EDAの3つの基本ステップ、(1) 生データの確認、(2) 要約統計量の算出、(3) データ可視化、に従って分析を進めます。

1. 生データの確認:

# データフレームの構造と型情報を表示
print("データ構造:")
tips.info()

# 先頭5行を表示
print("\nデータプレビュー (先頭):")
print(tips.head())

2. 要約統計量の算出:

# 数値変数の要約統計量
print("\n数値変数の要約統計量:")
print(tips.describe())

# カテゴリ変数の度数分布
print("\n性別ごとのサンプル数:")
print(tips['sex'].value_counts())

3. データの可視化: tip と total_bill の関係性を、sex ごとに色分けして可視化します。seabornライブラリの lmplot を使うと、散布図と回帰直線を簡潔に描画できます。

# seabornを用いて、性別で層別化した散布図と回帰直線をプロット
sns.lmplot(data=tips, x='total_bill', y='tip', hue='sex',
           height=6, aspect=1.5, palette={'Male': 'royalblue', 'Female': 'tomato'})

plt.title('チップ額と会計総額の関係（性別による層別化）')
plt.xlabel('会計総額 (Total Bill, $)')
plt.ylabel('チップ額 (Tip, $)')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

図6.1から、会計総額が増えるにつれてチップ額も増えるという、両性別で共通の正の傾向が見られます。しかし、2本の回帰直線の傾きはわずかに異なっているようにも見えます。この「傾きの違い」をモデルに組み込むのが交互作用モデルです。

6.1.2 交互作用モデル (Interaction Model)

図6.1に示された2本の異なる回帰直線を、1つの統一された回帰モデルで表現してみましょう。このモデルは、会計総額の効果が性別によって異なる可能性を捉えるため、交互作用項 (interaction term) を含みます。

まず、カテゴリ変数 sex を数値に変換（ダミー変数化）し、交互作用項 total_bill * is_male を作成します。

# ダミー変数を作成 (sex_Maleは男性なら1, 女性なら0)
# drop_first=Trueで、参照カテゴリ(Female)の列を削除
tips_encoded = pd.get_dummies(tips, columns=['sex'], drop_first=True)
tips_encoded.rename(columns={'sex_Male': 'is_male'}, inplace=True)

# 交互作用項 (total_bill * is_male) を作成
tips_encoded['bill_x_male'] = tips_encoded['total_bill'] * tips_encoded['is_male']

# 説明変数Xと目的変数yを定義
X = tips_encoded[['total_bill', 'is_male', 'bill_x_male']]
y = tips_encoded['tip']

# 交互作用モデルの学習
interaction_model = LinearRegression()
interaction_model.fit(X, y)

# 係数を表示
intercept = interaction_model.intercept_
coefs = interaction_model.coef_

print("交互作用モデルの係数:")
print(f"  切片 (b_0): {intercept:.4f}")
print(f"  total_bill (b_bill): {coefs[0]:.4f}")
print(f"  is_male (b_male): {coefs[1]:.4f}")
print(f"  bill_x_male (b_bill_male): {coefs[2]:.4f}")

このモデルは以下の数式で表現されます。

\[\widehat{\text{tip}} = b_0 + b_{\text{bill}} \cdot \text{total\_bill} + b_{\text{male}} \cdot \mathbb{1}_{\text{male}} + b_{\text{bill,male}} \cdot (\text{total\_bill} \times \mathbb{1}_{\text{male}})\]

ここで、\(\\mathbb{1}\_{\\text{male}}\) は指示関数 (indicator function) であり、客が男性の場合に1、女性の場合に0を取ります。

この式から、性別ごとの回帰直線を導出できます。

女性客の場合 (\(\mathbb{1}_{\text{male}} = 0\)):

\[ \begin{aligned} \widehat{\text{tip}}*{\text{female}} &= b_0 + b*{\text{bill}} \cdot \text{total\_bill} \\ &= 0.9203 + 0.1051 \cdot \text{total\_bill} \end{aligned} \]

\(b_0\) は女性客の切片、\(b_{\text{bill}}\) は女性客の会計総額に対する傾きを表します。
男性客の場合 (\(\mathbb{1}_{\text{male}} = 1\)):

\[ \begin{aligned} \widehat{\text{tip}}*{\text{male}} &= (b_0 + b*{\text{male}}) + (b_{\text{bill}} + b_{\text{bill,male}}) \cdot \text{total\_bill} \\ &= (0.9203 + 0.1349) + (0.1051 - 0.0003) \cdot \text{total\_bill} \\ &= 1.0552 + 0.1048 \cdot \text{total\_bill} \end{aligned} \]

男性客の切片は \((b_0 + b_{\text{male}})\)、傾きは \((b_{\text{bill}} + b_{\text{bill,male}})\) となります。\(b_{\text{male}}\) と \(b_{\text{bill,male}}\) は、それぞれ女性（参照カテゴリ）に対する切片と傾きの差分（オフセット）と解釈できます。

このモデルでは、会計総額がチップ額に与える影響（傾き）が性別によってわずかに異なりますが、ほぼ同じです。これが交互作用です。

6.1.3 平行傾斜モデル (Parallel Slopes Model)

もう一つのアプローチとして、会計総額の効果（傾き）は男女で共通であると仮定し、切片のみが異なるとするモデルを考えることができます。これを平行傾斜モデルと呼びます。このモデルは、交互作用項を含みません。

# 説明変数Xから交互作用項を除外
X_parallel = tips_encoded[['total_bill', 'is_male']]
y_parallel = tips_encoded['tip']

# 平行傾斜モデルの学習
parallel_model = LinearRegression()
parallel_model.fit(X_parallel, y_parallel)

# 係数を表示
intercept_p = parallel_model.intercept_
coefs_p = parallel_model.coef_

print("平行傾斜モデルの係数:")
print(f"  切片 (b_0): {intercept_p:.4f}")
print(f"  total_bill (b_bill): {coefs_p[0]:.4f}")
print(f"  is_male (b_male): {coefs_p[1]:.4f}")

モデルの数式は以下の通りです。 \(\(\widehat{\text{tip}} = b_0 + b_{\text{bill}} \cdot \text{total\_bill} + b_{\text{male}} \cdot \mathbb{1}_{\text{male}}\)\) 両者の傾きは共通で 0.1050 となっており、2本の直線は平行になります。

6.1.4 モデルの比較と解釈

交互作用モデルと平行傾斜モデルを視覚的に比較してみましょう。

# 交互作用モデルと平行傾斜モデルを並べてプロット
fig, axes = plt.subplots(1, 2, figsize=(16, 7), sharey=True)

# 交互作用モデルのプロット
sns.scatterplot(data=tips, x='total_bill', y='tip', hue='sex', alpha=0.6, ax=axes[0],
                palette={'Male': 'royalblue', 'Female': 'tomato'}, legend=False)
bill_range = np.array([tips['total_bill'].min(), tips['total_bill'].max()])
# 女性の直線
axes[0].plot(bill_range, intercept + coefs[0] * bill_range, color='tomato', linewidth=2)
# 男性の直線
axes[0].plot(bill_range, (intercept + coefs[1]) + (coefs[0] + coefs[2]) * bill_range, color='royalblue', linewidth=2)
axes[0].set_title('交互作用モデル', fontsize=14)
axes[0].set_xlabel('会計総額 ($)'); axes[0].set_ylabel('チップ額 ($)')
axes[0].grid(True, linestyle='--', alpha=0.6)

# 平行傾斜モデルのプロット
sns.scatterplot(data=tips, x='total_bill', y='tip', hue='sex', alpha=0.6, ax=axes[1],
                palette={'Male': 'royalblue', 'Female': 'tomato'})
# 女性の直線
axes[1].plot(bill_range, intercept_p + coefs_p[0] * bill_range, color='tomato', linewidth=2)
# 男性の直線
axes[1].plot(bill_range, (intercept_p + coefs_p[1]) + coefs_p[0] * bill_range, color='royalblue', linewidth=2)
axes[1].set_title('平行傾斜モデル', fontsize=14)
axes[1].set_xlabel('会計総額 ($)'); axes[1].set_ylabel('')
axes[1].grid(True, linestyle='--', alpha=0.6)
axes[1].legend(title='性別')

plt.suptitle('モデル比較: 交互作用モデル vs 平行傾斜モデル', fontsize=16)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()

このデータに関しては、両モデルの見た目の違いは非常に小さいです。傾きの違いがほとんどないため、交互作用モデルと平行傾斜モデルは酷似した結果を与えます。モデル選択については、6.3.1節でより詳しく議論します。

6.2 2つの数値変数を持つモデル

次に、2つの数値説明変数を持つ重回帰モデルを考えます。ここでは、seabornに付属するpenguinsデータセットを用います。

以下の回帰モデルを構築します。

目的変数 \(y\): ペンギンの体重 (body_mass_g)
説明変数:
1. \(x\_1\): くちばしの長さ (bill_length_mm)
2. \(x\_2\): フリッパーの長さ (flipper_length_mm)

6.2.1 探索的データ分析

データを読み込み、欠損値を持つ行を除外します。

# seabornからpenguinsデータセットをロード
penguins = sns.load_dataset("penguins")

# 欠損値を含む行を削除して分析対象を準備
penguins_ch6 = penguins[['body_mass_g', 'bill_length_mm', 'flipper_length_mm']].dropna()

# データプレビューと要約統計量
print(penguins_ch6.head())
print("\n", penguins_ch6.describe())

次に、3つの変数間の相関関係を相関行列 (correlation matrix) とヒートマップ (heatmap) で確認します。

# 相関行列を計算
correlation_matrix = penguins_ch6.corr()
print("\n相関行列:")
print(correlation_matrix)

# 相関行列をヒートマップで可視化
plt.figure(figsize=(7, 5))
sns.heatmap(correlation_matrix, annot=True, cmap='viridis', vmin=0, vmax=1, fmt=".3f")
plt.title('変数間の相関ヒートマップ')
plt.show()

ヒートマップから、体重 (body_mass_g) はフリッパーの長さ (flipper_length_mm) と非常に強い正の相関 (\(r=0.873\)) を、くちばしの長さ (bill_length_mm) とも中程度の正の相関 (\(r=0.595\)) を持つことがわかります。

6.2.2 回帰平面 (Regression Plane)

2つの数値説明変数を持つ回帰モデルは、3次元空間における回帰平面として表現されます。

# 重回帰モデルの構築
X = penguins_ch6[['bill_length_mm', 'flipper_length_mm']]
y = penguins_ch6['body_mass_g']

penguin_model = LinearRegression()
penguin_model.fit(X, y)

# 係数を表示
intercept_p = penguin_model.intercept_
coefs_p = penguin_model.coef_

print("回帰係数:")
print(f"  切片 (b_0): {intercept_p:.4f}")
print(f"  bill_length_mm (b_bill): {coefs_p[0]:.4f}")
print(f"  flipper_length_mm (b_flipper): {coefs_p[1]:.4f}")

得られた回帰平面の方程式は以下のようになります。 \(\(\widehat{\text{weight}} = -5980.6869 + 33.3499 \cdot \text{bill\_length} + 44.9101 \cdot \text{flipper\_length}\)\)

係数の解釈は以下の通りです。

b_bill (33.3499): フリッパーの長さを一定に保った場合、くちばしが1mm長くなると、体重は平均して約33.3g重くなると予測される。
b_flipper (44.9101): くちばしの長さを一定に保った場合、フリッパーが1mm長くなると、体重は平均して約44.9g重くなると予測される。

この回帰平面を3Dで可視化してみましょう。

# 3Dプロットの作成
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')

# 散布図をプロット
ax.scatter(penguins_ch6['bill_length_mm'], penguins_ch6['flipper_length_mm'], penguins_ch6['body_mass_g'],
           c='darkcyan', marker='o', alpha=0.5, label='観測データ')

# 回帰平面のグリッドを作成
bill_range = np.linspace(penguins_ch6['bill_length_mm'].min(), penguins_ch6['bill_length_mm'].max(), 30)
flipper_range = np.linspace(penguins_ch6['flipper_length_mm'].min(), penguins_ch6['flipper_length_mm'].max(), 30)
bill_grid, flipper_grid = np.meshgrid(bill_range, flipper_range)
mass_pred = intercept_p + coefs_p[0] * bill_grid + coefs_p[1] * flipper_grid

# 回帰平面をプロット
ax.plot_surface(bill_grid, flipper_grid, mass_pred, color='orange', alpha=0.4)

ax.set_xlabel('くちばしの長さ (mm)')
ax.set_ylabel('フリッパーの長さ (mm)')
ax.set_zlabel('体重 (g)')
ax.set_title('ペンギンの体重の回帰平面', fontsize=16)
ax.view_init(elev=20., azim=125)
plt.show()

6.3 関連トピック

6.3.1 モデル選択: 複雑さと適合度のトレードオフ

6.1節では、交互作用モデルと平行傾斜モデルの2つを比較しました。どちらのモデルを採用すべきでしょうか？このプロセスがモデル選択 (Model Selection) です。

指導原理はオッカムの剃刀 (Occam's Razor)、「他の条件が同じなら、より単純なモデルが望ましい」です。モデルの複雑さを増すのは、それがデータの構造を説明する上で実質的な改善をもたらす場合にのみ正当化されます。

tipsデータの例では、2つのモデルは酷似していました。これは交互作用の効果が非常に小さいためです。このような場合、交互作用項を含まない、より単純な平行傾斜モデルを選ぶのが合理的です。

6.3.2 決定係数 \(R^2\) によるモデル評価

モデルの適合度を定量的に評価する指標が決定係数 \(R^2\) (R-squared) です。\(R^2\) は、目的変数 \(y\) の全変動のうち、モデルによって説明された変動の割合を示します。

\[R^2 = \frac{Var(\hat{y})}{Var(y)} = 1 - \frac{Var(e)}{Var(y)}\]

\(R^2\) は 0 から 1 の範囲の値を取り、1に近いほどモデルがデータの変動を良く説明していることを意味します。

tipsデータで \(R^2\) を比較してみましょう。

# 交互作用モデルのR^2
r2_interaction = r2_score(y, interaction_model.predict(X))
# 平行傾斜モデルのR^2
r2_parallel = r2_score(y_parallel, parallel_model.predict(X_parallel))

print(f"交互作用モデル R²: {r2_interaction:.4f}")
print(f"平行傾斜モデル R²: {r2_parallel:.4f}")

\(R^2\) の差はごくわずかです。これは、交互作用項を追加してもモデルの説明力がほとんど向上しないことを示しており、より単純な平行傾斜モデルを選択するという判断を裏付けます。もし、交互作用モデルの\(R^2\)が大幅に高ければ、その複雑さは正当化され、交互作用モデルの採用が支持されます。

6.3.3 シンプソンのパラドックスと交絡変数

シンプソンのパラドックス (Simpson's Paradox) は、データ全体で観測される傾向が、データをサブグループに分割すると消滅、あるいは逆転する現象です。

この現象を penguins データで見てみましょう。ここでは、種 (species) が交絡変数 (confounding variable) の役割を果たす可能性があります。

まず、くちばしの長さ (bill_length_mm) と体重 (body_mass_g) の関係を全体で見てみます。次に、ペンギンの種（Adelie, Chinstrap, Gentoo）でグループ分けして関係を見ます。

# シンプソンのパラドックスを示すプロット
# species列を結合し、欠損値を除外
penguins_full = sns.load_dataset("penguins").dropna()

# lmplotで全体と層別化を同時に可視化
g = sns.lmplot(data=penguins_full, x='bill_length_mm', y='body_mass_g',
               hue='species', height=6, aspect=1.5, palette='viridis')

# 全体の回帰直線を追加
sns.regplot(data=penguins_full, x='bill_length_mm', y='body_mass_g',
            scatter=False, color='black', line_kws={'linestyle':'--'}, ax=g.ax)

g.ax.set_title('シンプソンのパラドックス: ペンギンの種による傾向の違い')
g.ax.set_xlabel('くちばしの長さ (mm)')
g.ax.set_ylabel('体重 (g)')
plt.legend(title='Species')
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

図6.5では、データ全体で見ると（黒い破線）、くちばしが長いほど体重が重いという明確な正の相関があります。しかし、種ごとに見ると（色のついた線）、特にアデリーペンギン（紫色）とジェンツーペンギン（黄色）のグループ内では、その傾きが全体よりも緩やかになっています。

これは、種によってくちばしの長さと体重の分布が大きく異なるためです。例えば、ジェンツーペンギンは他の種より全体的に体が大きく、くちばしも長いため、データ全体の関係性を「引き上げて」います。種という交絡変数を無視すると、くちばしの長さと体重の関係を過大評価してしまう可能性があります。

6.4 結論と次章への展望

本章では、seabornの標準データセットを用いて、重回帰分析の基礎を学びました。

重回帰: 複数の説明変数を用いて目的変数をモデル化する。
交互作用: ある説明変数の効果が、別の説明変数の水準によって変化する現象。
モデル選択: オッカムの剃刀を指針とし、モデルの複雑さと適合度（\(R^2\)など）のバランスを考える。
交絡: 第三の変数が原因で、2つの変数間に見せかけの相関が生じる現象（シンプソンのパラドックス）。

ここまでは、手元の標本 (sample) データに最もよくフィットする線を「推定」することに焦点を当ててきました。次章から始まるパートIII「統計的推測」では、この標本から得られた知見を、より大きな母集団 (population) 全体に一般化するための理論と手法を学びます。