Bài thực hành từ ISLP (James, Witten, Hastie, Tibshirani, 2023) — §3.6, tr. 123–128
Trong bài lab này, chúng ta sẽ áp dụng mọi thứ đã học ở Chương 3 vào dữ liệu thực tế bằng Python. Bạn sẽ dùng statsmodels và scikit-learn để ước lượng hồi quy tuyến tính, thêm biến tương tác, biến đổi phi tuyến, xử lý biến định tính, và tự viết hàm đánh giá mô hình.
在本實驗中，我們將把第 3 章所學的一切應用到真實數據上。你將使用 statsmodels 和 scikit-learn 來估計線性迴歸、添加交互項、非線性變換、處理定性變數，並自行編寫模型評估函數。
Chúng ta dùng bộ dữ liệu Boston (giá nhà) từ gói ISLP và Auto (dữ liệu xe hơi).
使用來自 ISLP 包的 Boston（房價）和 Auto（汽車數據）數據集。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import summarize, ModelSpec as MS
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# Tải dữ liệu
Boston = load_data('Boston')
Auto = load_data('Auto')Dùng statsmodels để ước lượng mối quan hệ giữa medv (giá nhà trung vị) và lstat (tỷ lệ dân thu nhập thấp).
使用 statsmodels 估計 medv（房價中位數）與 lstat（低收入人口比例）之間的關係。
# Cách 1: Dùng statsmodels (OLS chuẩn)
X = sm.add_constant(Boston['lstat']) # Thêm intercept
y = Boston['medv']
model = sm.OLS(y, X).fit()
print(model.summary().tables[1])
# Cách 2: Dùng sklearn
X_sk = Boston[['lstat']]
y_sk = Boston['medv']
lr = LinearRegression().fit(X_sk, y_sk)
print(f"β₀ = {lr.intercept_:.3f}, β₁ = {lr.coef_[0]:.3f}")Kết quả: β₁ ≈ −0.95, nghĩa là khi tỷ lệ dân thu nhập thấp tăng 1%, giá nhà giảm khoảng $950. p-value ≈ 0 → có ý nghĩa thống kê mạnh.
結果：β₁ ≈ −0.95，表示低收入人口比例每增加 1%，房價約下降 $950。p-value ≈ 0 → 統計上極顯著。
Thêm nhiều biến dự báo: lstat + age (tuổi nhà).
添加多個預測變數：lstat + age（房屋年齡）。
# Dùng tất cả biến
predictors = ['lstat', 'age']
X_multi = sm.add_constant(Boston[predictors])
model_multi = sm.OLS(y, X_multi).fit()
print(model_multi.summary().tables[1])
# Hoặc dùng tất cả 13 biến
X_all = sm.add_constant(Boston.drop('medv', axis=1))
model_all = sm.OLS(y, X_all).fit()
print(f"R² = {model_all.rsquared:.4f}, Adj R² = {model_all.rsquared_adj:.4f}")Lưu ý: R² luôn tăng khi thêm biến, nhưng R² hiệu chỉnh (Adjusted R²) sẽ giảm nếu biến thêm vào không có ích. Dùng Adjusted R² để so sánh mô hình.
注意：R² 添加變數時總是增加，但 調整後 R² 若新增變數無用則會下降。用調整後 R² 來比較模型。
Tạo biến tương tác giữa lstat và age bằng ModelSpec của ISLP.
使用 ISLP 的 ModelSpec 創建 lstat 和 age 之間的交互項。
# Dùng ModelSpec để tạo tương tác
design = MS(['lstat', 'age', ('lstat', 'age')]).fit(
Boston[['lstat', 'age']]
)
X_inter = design.transform(Boston[['lstat', 'age']])
model_inter = sm.OLS(y, X_inter).fit()
print(model_inter.summary().tables[1])Nếu hệ số tương tác có ý nghĩa thống kê (p < 0.05), điều đó cho thấy tác động của lstat lên medv phụ thuộc vào age.
若交互項係數具統計顯著性（p < 0.05），表示 lstat 對 medv 的影響取決於 age。
Thêm lstat² vào mô hình để nắm bắt quan hệ phi tuyến.
加入 lstat² 來捕捉非線性關係。
# Tạo biến bậc hai
Boston['lstat2'] = Boston['lstat'] ** 2
# Dùng ModelSpec với hàm tự định nghĩa
design_poly = MS([MS('lstat').poly(2), 'age']).fit(Boston)
X_poly = design_poly.transform(Boston)
model_poly = sm.OLS(y, X_poly).fit()
print(model_poly.summary().tables[1])
# Vẽ đường hồi quy
X_sort = Boston['lstat'].sort_values()
X_pred = sm.add_constant(
pd.DataFrame({'lstat': X_sort, 'lstat2': X_sort**2})
)
y_pred = model_poly.predict(X_pred)
plt.scatter(Boston['lstat'], y, alpha=0.5)
plt.plot(X_sort, y_pred, 'r-', linewidth=2)
plt.xlabel('lstat'); plt.ylabel('medv')MS(...).poly(degree) để tự động tạo các số hạng đa thức, thay vì tự tính thủ công.MS(...).poly(degree) 自動生成多項式項，而非手動計算。Dùng dữ liệu Carseats với biến ShelveLoc (vị trí kệ: Bad/Medium/Good).
使用 Carseats 數據，含 ShelveLoc（貨架位置：Bad/Medium/Good）。
Carseats = load_data('Carseats')
# ModelSpec tự động tạo dummy variables
design_cat = MS(['Price', 'Income', 'ShelveLoc']).fit(Carseats)
X_cat = design_cat.transform(Carseats)
y_cat = Carseats['Sales']
model_cat = sm.OLS(y_cat, X_cat).fit()
print(model_cat.summary().tables[1])
# Kết quả: ShelveLoc[Good] có hệ số dương lớn nhất →
# vị trí kệ tốt làm tăng doanh sốGiải thích: ModelSpec tự động tạo biến giả (dummy) cho biến định tính, bỏ qua một hạng mục làm baseline. Ở đây ShelveLoc[Bad] là baseline.
解釋：ModelSpec 自動為定性變數創建虛擬變數，省略一個類別作為基準。此處 ShelveLoc[Bad] 為基準。
Tạo hàm eval_mse() để tính MSE trên tập kiểm tra — điều mà statsmodels không làm sẵn.
創建 eval_mse() 函數來計算測試集 MSE——statsmodels 沒有內建此功能。
def eval_mse(model, X_test, y_test):
"""Tính MSE trên tập kiểm tra cho bất kỳ mô hình nào có .predict()"""
y_pred = model.predict(X_test)
return np.mean((y_test - y_pred) ** 2)
# Chia dữ liệu
X = Boston.drop('medv', axis=1)
y = Boston['medv']
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3, random_state=42
)
# Huấn luyện + đánh giá
model = LinearRegression().fit(X_train, y_train)
mse = eval_mse(model, X_test, y_test)
print(f"Test MSE = {mse:.2f}")| Kỹ thuật | Công cụ | Khi nào dùng |
|---|---|---|
| Hồi quy tuyến tính đơn | sm.OLS / LinearRegression | 1 biến dự báo, quan hệ tuyến tính |
| Hồi quy tuyến tính bội | sm.OLS + nhiều biến | Nhiều biến, cần diễn giải |
| Biến tương tác | MS(...) với tuple | Khi tác động của biến này phụ thuộc biến kia |
| Đa thức | MS(...).poly(degree) | Quan hệ phi tuyến nhẹ |
| Biến định tính | MS(...) tự động dummy | Biến phân loại (giới tính, vùng...) |
| Đánh giá | train_test_split + MSE | Luôn luôn! |
🔥 Bạn vừa hoàn thành Chương 3! Bạn đã có đủ công cụ để phân tích hồi quy trên dữ liệu thực tế. Ở Chương 4, chúng ta sẽ chuyển sang bài toán phân loại (classification).
| 技術 | 工具 | 使用時機 |
|---|---|---|
| 簡單線性迴歸 | sm.OLS / LinearRegression | 1 個預測變數，線性關係 |
| 多元線性迴歸 | sm.OLS + 多變數 | 多變數，需要可解釋性 |
| 交互項 | MS(...) 帶 tuple | 當某變數的影響取決於另一變數 |
| 多項式 | MS(...).poly(degree) | 輕微非線性關係 |
| 定性變數 | MS(...) 自動虛擬化 | 類別變數（性別、地區…） |
| 評估 | train_test_split + MSE | 永遠需要！ |
🔥 你剛完成了第 3 章！你已具備在真實數據上進行迴歸分析的所有工具。在第 4 章，我們將轉向分類問題。