8. Inference for Regression

ブートストラップ法による回帰係数の信頼区間推定

スライド

最終課題

1. 線形回帰モデルについて

まず、基本となる線形回帰モデルから始めましょう。線形回帰は、一つの変数（説明変数、\(x\)）を使って、もう一つの変数（目的変数、\(y\)）を予測・説明するための最も基本的なモデルです。

モデルは以下のような数式で表されます。

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

各項の意味は以下の通りです。

\(y_i\): i番目のサンプルの目的変数（例：家の価格）
\(x_i\): i番目のサンプルの説明変数（例：家の広さ）
\(\beta_0\): 切片 (Intercept)。\(x\)が0のときの\(y\)の予測値。
\(\beta_1\): 回帰係数 (Regression Coefficient) または 傾き (Slope)。\(x\)が1単位増加したときの\(y\)の予測値の変化量。私たちが最も関心を持つことが多いパラメータです。
\(\varepsilon_i\): 誤差項 (Error Term)。モデルでは説明しきれない部分（ノイズなど）。

私たちの目的は、手元にあるデータ \((x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\) を使って、真の（しかし未知の）パラメータである \(\beta_0\) と \(\beta_1\) を推定することです。この推定された値を、ハット記号を付けて \(\hat{\beta}_0, \hat{\beta}_1\) と書きます。

最も一般的な推定方法は最小二乗法 (Ordinary Least Squares, OLS) で、予測値と実際の値の差（残差）の二乗和を最小化することで、最適な \(\hat{\beta}_0\) と \(\hat{\beta}_1\) を見つけます。

Pythonコード例：基本的な線形回帰

まずは、綺麗な線形関係を持つデータを生成し、statsmodels を使って回帰モデルを当てはめてみましょう。

# 必要なライブラリをインポート
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns

# スタイルの設定
sns.set_style('whitegrid')

# 再現性のための乱数シード固定
np.random.seed(42)

# 1. データの生成
# 説明変数 X (0から10の範囲で100個)
X = np.linspace(0, 10, 100)
# 真の係数
beta_0_true = 2.0
beta_1_true = 5.0
# 正規分布に従う誤差
epsilon = np.random.normal(0, 2.0, 100)
# 目的変数 y
y = beta_0_true + beta_1_true * X + epsilon

# 2. statsmodelsのためのデータ準備
# statsmodelsでは切片項をモデルに含めるために定数項(const)を追加する必要がある
X_const = sm.add_constant(X)

# 3. モデルの学習 (最小二乗法)
model = sm.OLS(y, X_const)
results = model.fit()

# 4. 結果の表示
print(results.summary())

# 5. 結果の可視化
plt.figure(figsize=(10, 6))
plt.scatter(X, y, label='Data points')
plt.plot(X, results.predict(X_const), color='red', linewidth=2, label='OLS Regression Line')
plt.title('Basic Linear Regression')
plt.xlabel('X (Explanatory Variable)')
plt.ylabel('y (Dependent Variable)')
plt.legend()
plt.show()

results.summary() の出力を見ると、const ( \(\hat{\beta}_0\)) と x1 ( \(\hat{\beta}_1\)) の推定値 (coef) が、私たちが設定した真の値 2.0 と 5.0 に非常に近いことがわかります。

2. 線形回帰モデルの回帰係数のばらつきについて

results.summary() の結果には、coef (係数の推定値) だけでなく、std err (標準誤差) や [0.025, 0.975] (95%信頼区間) といった列も表示されています。これらはなぜ必要なのでしょうか？

それは、私たちが得た回帰係数 \(\hat{\beta}_1\) は、あくまで手元にある「サンプルデータ」から計算された推定値に過ぎないからです。もし、全く別のサンプルデータを母集団から取ってくれば、計算される \(\hat{\beta}_1\) は少し違う値になるでしょう。

この「もし別のサンプルだったら、係数はどれくらい変動するだろうか？」というばらつきの大きさを表すのが標準誤差 (Standard Error) であり、そのばらつきを考慮して「真の \(\beta_1\) が含まれるであろう区間」を推定したものが信頼区間 (Confidence Interval) です。

statsmodels が計算する標準誤差や信頼区間は、誤差項 \(\varepsilon_i\) が以下の理論的な仮定を満たすことを前提としています。

正規性: 誤差は正規分布に従う。
等分散性: 誤差の分散は、\(x\) の値によらず一定である。
独立性: 誤差は互いに独立である。

しかし、現実のデータではこれらの仮定が満たされないことも少なくありません。特に、次に示す外れ値が存在する場合、係数の推定は不安定になります。

3. 線形回帰モデルの外れ値に対する回帰係数の不安定性

最小二乗法（OLS）は、残差の「二乗」和を最小化するため、一つでも非常に大きな残差（つまり、回帰直線から大きく外れた点）があると、その点に強く影響を受けてしまいます。

先ほどのデータに意図的に外れ値を一つ加えて、回帰係数がどう変わるか見てみましょう。

Pythonコード例：外れ値の影響

# 1. 外れ値の追加
# 元のデータに外れ値を追加
X_outlier = np.append(X, 1)   # x=1 の位置に
y_outlier = np.append(y, 40)  # y=40 という大きな値を追加

# 2. statsmodelsのためのデータ準備
X_outlier_const = sm.add_constant(X_outlier)

# 3. モデルの再学習
model_outlier = sm.OLS(y_outlier, X_outlier_const)
results_outlier = model_outlier.fit()

# 4. 結果の表示
print("--- Original Model ---")
print(f"Coefficient (beta_1): {results.params[1]:.4f}")
print(f"95% CI: [{results.conf_int().loc[1][0]:.4f}, {results.conf_int().loc[1][1]:.4f}]")
print("\n--- Model with Outlier ---")
print(results_outlier.summary())


# 5. 可視化による比較
plt.figure(figsize=(12, 7))
# 元のデータと回帰直線
plt.scatter(X, y, label='Original Data')
plt.plot(X, results.predict(X_const), color='blue', linestyle='--', label=f'Original Line (beta_1={results.params[1]:.2f})')
# 外れ値を含むデータと回帰直線
plt.scatter(X_outlier, y_outlier, label='Data with Outlier', alpha=0.6)
plt.plot(X_outlier, results_outlier.predict(X_outlier_const), color='red', linewidth=2, label=f'Outlier Line (beta_1={results_outlier.params[1]:.2f})')

plt.title('Impact of an Outlier on Linear Regression')
plt.xlabel('X')
plt.ylabel('y')
plt.legend()
plt.show()

結果の考察: たった一つの外れ値を加えただけで、回帰係数 \(\hat{\beta}_1\) は 5.04 から 3.66 へと劇的に変化しました。回帰直線が外れ値に「引っ張られて」いるのが視覚的にも明らかです。

このように、OLSの理論的な信頼区間は、外れ値のような仮定を乱すデータがあると、信頼性が低くなります。そこで登場するのがブートストラップ法です。ブートストラップは、難しい数理モデルの仮定に頼らず、データ自身の力で係数のばらつきを推定する非常に強力な手法です。

4. ブートストラップ法の解説

ブートストラップ法は、手元にあるサンプルデータを「ミニチュアの母集団」と見なして、そこから何度もシミュレーション（リサンプリング）を行うことで、統計量のばらつきを経験的に調べる手法です。

アルゴリズムは驚くほどシンプルです。

リサンプリング: 元のデータセット（サイズ \(N\)）から、重複を許して (with replacement) \(N\) 個のデータをランダムに抽出します。こうして作られた新しいデータセットを「ブートストラップサンプル」と呼びます。
- 重複を許すため、元のデータのある点は複数回選ばれるかもしれず、一度も選ばれない点も出てきます。
統計量の計算: 作成したブートストラップサンプルに対して、興味のある統計量（今回は回帰係数 \(\hat{\beta}_1\)）を計算します。
繰り返し: 上記のステップ1と2を、非常に多数回（例えば、B=1000回や10000回）繰り返します。

この結果、私たちは B個の \(\hat{\beta}_1\) の推定値の集まり（分布）を得ることができます。この分布が、真のサンプリング分布の近似となるのです。この分布のばらつきを見ることで、元の係数推定値の不確実性を評価できます。

ブートストラップ法の最大の利点は、正規分布などの強い仮定を必要としない点です。データが持つ本来の分布構造（外れ値を含むなど）をリサンプリングによって再現するため、より現実に即した頑健な（ロバストな）推定が可能になります。

5. 線形回帰モデルに対する回帰係数のブートストラップ法による区間推定

それでは、先ほどの外れ値を含むデータセットに対して、ブートストラップ法を適用し、回帰係数 \(\hat{\beta}_1\) の信頼区間を求めてみましょう。

Pythonコード例：ブートストラップ法の実装

# 外れ値を含むデータセットを使用
# X_outlier, y_outlier

# 1. ブートストラップの設定
n_iterations = 1000  # ブートストラップの繰り返し回数 B
n_size = len(X_outlier) # サンプルサイズ N

# 2. 係数を保存するためのリストを初期化
bootstrap_betas = []

# 3. ブートストラップ・ループ
for i in range(n_iterations):
    # a. リサンプリング (重複を許してインデックスを抽出)
    indices = np.random.choice(range(n_size), size=n_size, replace=True)

    # b. ブートストラップサンプルを作成
    X_sample = X_outlier[indices]
    y_sample = y_outlier[indices]

    # c. ブートストラップサンプルに対してモデルを学習
    X_sample_const = sm.add_constant(X_sample)
    model_boot = sm.OLS(y_sample, X_sample_const)
    results_boot = model_boot.fit()

    # d. 回帰係数(beta_1)を保存
    bootstrap_betas.append(results_boot.params[1])

# 4. 結果の確認
# リストをnumpy配列に変換
bootstrap_betas = np.array(bootstrap_betas)

# ブートストラップで得られた係数の分布を可視化
plt.figure(figsize=(10, 6))
sns.histplot(bootstrap_betas, kde=True)
plt.title('Distribution of Bootstrap Coefficients (beta_1)')
plt.xlabel('Estimated Coefficient value')
plt.ylabel('Frequency')
plt.show()

# 5. パーセンタイル法による95%信頼区間の計算
# 2.5パーセンタイルと97.5パーセンタイルを求める
alpha = 0.05
lower_bound = np.percentile(bootstrap_betas, 100 * (alpha / 2))
upper_bound = np.percentile(bootstrap_betas, 100 * (1 - alpha / 2))

print(f"Bootstrap 95% Confidence Interval for beta_1:")
print(f"[{lower_bound:.4f}, {upper_bound:.4f}]")

# OLSの理論的な信頼区間と比較
print("\nFor comparison:")
print(f"OLS Theoretical 95% CI (with outlier):")
ci_outlier = results_outlier.conf_int().loc[1]
print(f"[{ci_outlier[0]:.4f}, {ci_outlier[1]:.4f}]")

6. ブートストラップによって得られた区間の解釈

上記のコードを実行すると、ブートストラップによって得られた多数の係数（今回は1000個）の分布がヒストグラムで表示されます。この分布が、データ自身から推定した「回帰係数のありえそうな値の範囲」です。

そして、この分布の下側2.5%点と上側2.5%点を結んで計算された区間がブートストラップ95%信頼区間です。

解釈: 「もし、私たちがこの実験（データ収集と分析）を何度も繰り返すことができたなら、そのうちの95%のケースで、構築されるこの信頼区間が『真の』回帰係数 \(\beta_1\) を含むだろう」と解釈します。

これは、「この区間内に95%の確率で真の値がある」という単純な確率論とは少し異なりますが、実用上は「真の係数がこの範囲にあると95%の信頼度で推定される」と考えて差し支えありません。

OLSの理論的信頼区間との比較: 外れ値があるデータでは、多くの場合、ブートストラップ信頼区間はOLSの理論的信頼区間よりも広くなります。これは、ブートストラップがリサンプリングの過程で外れ値の影響による係数の「ブレ」を正直に反映しているためです。これにより、私たちは係数推定の不確実性について、より現実的で保守的な評価を得ることができます。

まとめ

本講義では、以下の内容を学びました。

線形回帰の係数はサンプルデータに依存するため、常に**不確実性（ばらつき）**を伴います。
OLS（最小二乗法）は外れ値に非常に敏感で、係数の推定値が大きく歪められる可能性があります。
ブートストラップ法は、難しい理論的仮定に頼らず、データからのリサンプリングによって係数のばらつきを経験的に推定する強力な手法です。
ブートストラップ信頼区間は、特にデータに外れ値があったり、誤差の正規性が疑わしい場合に、より頑健な（ロバストな）推定を与えてくれます。

データ分析を行う際には、単にモデルを当てはめて係数を報告するだけでなく、その係数がどれほど信頼できるのかを評価することが極めて重要です。ブートストラップ法は、そのための強力なツールの一つとして、ぜひ覚えておきましょう。