3.2 多元線性迴歸

從簡單迴歸到多元迴歸：為什麼需要多個預測變數？

在 §3.1 中，我們用「電視廣告預算」一個變數來預測銷售額，得到了不錯的結果。但現實世界遠比這複雜——銷售額可能同時受到電視、收音機、報紙等多個廣告管道的影響。如果我們分別對每個預測變數擬合一個簡單線性迴歸，會產生兩個嚴重問題：

解決方案是多元線性迴歸模型（Multiple Linear Regression）——將所有預測變數放入同一個模型，每個變數擁有自己的斜率係數。

其中 \(\beta_j\) 代表在其他預測變數保持不變的情況下，\(X_j\) 每增加一單位對 \(Y\) 的平均影響。以廣告資料為例：

3.2.1 估計迴歸係數：最小平方法的推廣

與簡單線性迴歸類似，我們使用最小平方法（Least Squares）來估計 \(\beta_0, \beta_1, \ldots, \beta_p\)——選擇能使殘差平方和（RSS）最小化的參數值：

不同於簡單線性迴歸有封閉解（closed-form），多元迴歸的係數估計需要用矩陣代數來表示，形式較為複雜，但任何統計軟體都能輕鬆計算。在幾何上，當 \(p=2\) 時，最小平方解對應到三維空間中的一個平面（而非簡單迴歸的直線），該平面最小化了所有觀測點到平面的垂直距離平方和。

廣告資料的實證結果

使用 TV、radio、newspaper 三個變數預測 sales，得到以下係數估計（表 3.4）：

3.2.2 多元迴歸的四大核心問題

進行多元線性迴歸時，我們通常關心以下四個問題：

1. 至少有一個預測變數對反應變數有預測力嗎？
2. 所有預測變數都有幫助，還是只需要部分子集？
3. 模型對資料的擬合程度如何？
4. 給定一組預測變數值，該如何預測反應變數？預測有多準確？

問題一：反應變數與預測變數之間是否存在關係？

在簡單迴歸中，檢查 \(\beta_1 = 0\) 即可。在多元迴歸中，我們需要同時檢驗所有係數是否皆為零：

為此我們使用 F 統計量（F-statistic）：

其中 TSS（Total Sum of Squares）= \(\sum (y_i - \bar{y})^2\)。如果線性模型假設成立且 \(H_0\) 為真，則 F 統計量的期望值接近 1；若 \(H_a\) 為真，F 會顯著大於 1。

此外，F 檢定也可以檢驗特定子集的係數是否為零。假設我們想檢驗最後 q 個變數是否需要（即 \(H_0: \beta_{p-q+1} = \cdots = \beta_p = 0\)），可以擬合不包含這 q 個變數的簡化模型（其殘差平方和為 \(\text{RSS}_0\)），然後計算：

問題二：哪些變數真正重要？——變數選擇

當 F 檢定顯示「至少有一個預測變數有用」後，下一步是找出哪些變數才是真正重要的。這就是變數選擇（Variable Selection）問題（將在第六章深入探討）。

三種經典的自動化變數選擇策略：

模型品質的評估指標包括：Mallow's \(C_p\)、AIC（Akaike Information Criterion）、BIC（Bayesian Information Criterion）、以及調整後 \(R^2\)（Adjusted \(R^2\)）。這些指標在擬合度和模型複雜度之間取得平衡。

問題三：模型擬合度——RSE 與 R²

在多元迴歸中，我們同樣使用 RSE 和 \(R^2\) 來評估模型擬合度，但計算方式需要調整：

由於 \(R^2\) 在加入更多變數時永遠不會下降（即使加入的是雜訊變數），單純用 \(R^2\) 來比較不同大小的模型會傾向選擇最複雜的模型。這就是為什麼在變數選擇中，我們偏好使用調整後 \(R^2\)、AIC、BIC 等懲罰模型複雜度的指標。

問題四：預測與不確定性

一旦擬合了多元迴歸模型，就可以用 (3.21) 來做預測：

但預測伴隨著三種不確定性來源：

Python 實作：多元線性迴歸完整示範

# 3.2 多元線性迴歸 — 完整可執行示範
# 資料擷取與模型擬合
try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
from scipy import stats

# ===== 1. 載入 Advertising 資料 =====
import urllib.request, os

csv_path = os.path.join(DATA_PATH, 'Advertising.csv')
if not os.path.exists(csv_path):
    os.makedirs(DATA_PATH, exist_ok=True)
    url = 'https://www.statlearning.com/s/Advertising.csv'
    urllib.request.urlretrieve(url, csv_path)

ad = pd.read_csv(csv_path)
print("資料維度:", ad.shape)
print(ad.head())

# ===== 2. 擬合多元線性迴歸 (TV + radio + newspaper) =====
X = ad[['TV', 'radio', 'newspaper']].values
y = ad['sales'].values

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

print("\n===== 多元迴歸係數 =====")
print(f"Intercept:  {model.intercept_:.4f}")
print(f"TV:         {model.coef_[0]:.4f}")
print(f"radio:      {model.coef_[1]:.4f}")
print(f"newspaper:  {model.coef_[2]:.4f}")

# ===== 3. 模型評估 =====
y_pred = model.predict(X)
n, p = X.shape

rss = np.sum((y - y_pred) ** 2)
tss = np.sum((y - y.mean()) ** 2)
rse = np.sqrt(rss / (n - p - 1))
r2 = 1 - rss / tss
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)

# F-statistic
f_stat = ((tss - rss) / p) / (rss / (n - p - 1))
f_pvalue = 1 - stats.f.cdf(f_stat, p, n - p - 1)

print(f"\n===== 模型擬合統計 =====")
print(f"RSS:          {rss:.4f}")
print(f"RSE:          {rse:.4f}")
print(f"R²:           {r2:.4f}")
print(f"Adjusted R²:  {adj_r2:.4f}")
print(f"F-statistic:  {f_stat:.2f}")
print(f"F p-value:    {f_pvalue:.2e}")

# ===== 4. 比較：只用 TV 的簡單迴歸 vs. 多元迴歸 =====
model_tv = LinearRegression()
model_tv.fit(ad[['TV']].values, y)
y_pred_tv = model_tv.predict(ad[['TV']].values)
rss_tv = np.sum((y - y_pred_tv) ** 2)
r2_tv = r2_score(y, y_pred_tv)

model_tv_radio = LinearRegression()
model_tv_radio.fit(ad[['TV', 'radio']].values, y)
y_pred_tr = model_tv_radio.predict(ad[['TV', 'radio']].values)
rss_tr = np.sum((y - y_pred_tr) ** 2)
r2_tr = r2_score(y, y_pred_tr)

print(f"\n===== 模型比較 =====")
print(f"{'模型':<25} {'RSS':>10} {'R²':>8}")
print(f"{'TV only':<25} {rss_tv:>10.2f} {r2_tv:>8.4f}")
print(f"{'TV + radio':<25} {rss_tr:>10.2f} {r2_tr:>8.4f}")
print(f"{'TV + radio + newspaper':<25} {rss:>10.2f} {r2:>8.4f}")

# ===== 5. 預測：信賴區間 vs. 預測區間 =====
# 使用 statsmodels 來計算區間（sklearn 不直接支援）
import statsmodels.api as sm
X_sm = sm.add_constant(X)
model_sm = sm.OLS(y, X_sm).fit()

# 新資料點：TV=100, radio=20, newspaper=30
X_new = np.array([[1, 100, 20, 30]])
pred = model_sm.get_prediction(X_new)
pred_summary = pred.summary_frame(alpha=0.05)

print(f"\n===== 預測 (TV=100, radio=20, newspaper=30) =====")
print(f"點估計:           {pred_summary['mean'][0]:.2f}")
print(f"95% 信賴區間:     [{pred_summary['obs_ci_lower'][0]:.2f}, {pred_summary['obs_ci_upper'][0]:.2f}]")
print(f"95% 預測區間:     [{pred_summary['mean_ci_lower'][0]:.2f}, {pred_summary['mean_ci_upper'][0]:.2f}]")
# Note: statsmodels uses 'obs_ci' for prediction intervals and 'mean_ci' for confidence intervals

# ===== 6. 視覺化：實際 vs. 預測 =====
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# 6a. 實際 vs 預測散點圖
axes[0].scatter(y, y_pred, alpha=0.6, edgecolors='k', linewidth=0.5)
axes[0].plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2, label='完美預測')
axes[0].set_xlabel('實際銷售額')
axes[0].set_ylabel('預測銷售額')
axes[0].set_title(f'實際 vs. 預測 (R²={r2:.3f})')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 6b. 殘差圖
residuals = y - y_pred
axes[1].scatter(y_pred, residuals, alpha=0.6, edgecolors='k', linewidth=0.5)
axes[1].axhline(y=0, color='r', linestyle='--', lw=2)
axes[1].set_xlabel('預測銷售額')
axes[1].set_ylabel('殘差')
axes[1].set_title('殘差圖')
axes[1].grid(True, alpha=0.3)

# 6c. 係數比較圖
coef_names = ['TV', 'radio', 'newspaper']
coef_vals = model.coef_
colors_bar = ['#58a6ff' if v > 0 else '#f85149' for v in coef_vals]
axes[2].barh(coef_names, coef_vals, color=colors_bar, edgecolor='white')
axes[2].axvline(x=0, color='white', lw=1)
axes[2].set_xlabel('係數值')
axes[2].set_title('迴歸係數比較')
for i, v in enumerate(coef_vals):
    axes[2].text(v + 0.002, i, f'{v:.4f}', va='center', fontsize=9)

plt.tight_layout()
plt.savefig(os.path.join(DATA_PATH, '3_2_multiple_regression.png'), dpi=120, bbox_inches='tight')
plt.show()
print(f"\n✅ 圖片已儲存至: {os.path.join(DATA_PATH, '3_2_multiple_regression.png')}")

方法比較：簡單迴歸 vs. 多元迴歸