Bài 19: Multiple Linear Regression — nhiều feature

Bài 15 đã giới thiệu Linear Regression đơn biến — 1 feature dự đoán 1 target:

\[ \hat{y} = w x + b \]

Trong đó w là slope và b là intercept. Bài 16 cài bằng sklearn, Bài 17 đo bằng MSE/RMSE, Bài 18 đánh giá bằng R². Bài này mở rộng lên trường hợp thực tế: dataset có nhiều feature.

Vì sao 1 feature không đủ? Hiếm có hiện tượng nào trong thực tế chỉ phụ thuộc một biến. Giá nhà phụ thuộc cùng lúc vào diện tích, số phòng, vị trí, tuổi nhà. Lương phụ thuộc kinh nghiệm, ngành, công ty, thành phố. Bỏ feature liên quan ra ngoài → model underfit và coefficient của feature còn lại bị thiên lệch (omitted variable bias).

Với n feature \( x_1, x_2, \ldots, x_n \), mỗi feature có weight riêng \( w_1, w_2, \ldots, w_n \) và toàn model vẫn chỉ có 1 bias b:

\[ \hat{y} = w_1 x_1 + w_2 x_2 + \ldots + w_n x_n + b \]

Đặt vector weight \( \mathbf{w} = [w_1, w_2, \ldots, w_n]^T \) và vector feature \( \mathbf{x} = [x_1, x_2, \ldots, x_n]^T \). Phương trình thu gọn thành dot product:

\[ \hat{y} = \mathbf{w}^T \mathbf{x} + b \]

Đây vẫn là cùng một công thức của Bài 15, chỉ khác ở chỗ \( \mathbf{w}, \mathbf{x} \) bây giờ là vector chiều n thay vì scalar. Tính chất tuyến tính theo parameter giữ nguyên — model vẫn là Linear Regression.

Lưu ý ký hiệu: \( \mathbf{w}^T \mathbf{x} \) đôi khi viết là \( \mathbf{w} \cdot \mathbf{x} \) (dot product), trong code NumPy là w @ x. Ba cách viết tương đương.

Công thức trên áp dụng cho 1 sample. Trong thực tế bạn predict đồng thời cho cả batch n sample (lưu ý: cùng ký hiệu n với ngữ cảnh khác — ở đây n là số sample, d là số feature).

Gom batch lại thành ma trận thiết kế \( X \) shape \( (n, d) \): mỗi hàng là một sample, mỗi cột là một feature. Target gom thành vector \( \mathbf{y} \) shape \( (n,) \). Prediction cho cả batch:

\[ \hat{\mathbf{y}} = X \mathbf{w} + b \]

Với:

• \( X \) shape \( (n, d) \) — feature matrix.
• \( \mathbf{w} \) shape \( (d,) \) — weight vector.
• \( b \) scalar — bias (broadcast cộng vào mọi sample).
• \( \hat{\mathbf{y}} \) shape \( (n,) \) — vector prediction.

Mẹo gộp bias vào weight: thêm cột toàn 1 vào X (shape \( (n, d+1) \)) và đặt \( \mathbf{w}' = [w_1, \ldots, w_d, b]^T \). Lúc này:

\[ \hat{\mathbf{y}} = X \mathbf{w}' \]

Không còn bias riêng — quy ước này tiện cho dẫn xuất OLS ở mục 4.

Loss MSE viết dưới dạng ma trận (sau khi đã gộp bias vào X):

\[ L(\mathbf{w}) = \frac{1}{n} \| \mathbf{y} - X \mathbf{w} \|^2 \]

Lấy gradient theo \( \mathbf{w} \), đặt bằng 0:

\[ \nabla_{\mathbf{w}} L = -\frac{2}{n} X^T (\mathbf{y} - X \mathbf{w}) = 0 \]

Giải ra normal equation:

\[ \hat{\mathbf{w}} = (X^T X)^{-1} X^T \mathbf{y} \]

Công thức này đã preview ở Series 1 B27. Nó cho nghiệm OLS đúng, không cần iterate. Điều kiện: \( X^T X \) phải khả nghịch — đồng nghĩa các cột của X độc lập tuyến tính.

Khi nào ma trận singular hoặc gần singular:

• Hai feature trùng nhau hoàn toàn (vd "diện tích m²" và "diện tích ft²" cùng đơn vị quy đổi).
• One-Hot Encoding không drop cột nào — cột bias (toàn 1) bằng tổng các cột OHE.
• Số sample \( n \) ít hơn số feature \( d \) — \( X^T X \) chắc chắn singular.

Sklearn tránh bài toán nghịch đảo trực tiếp bằng SVD — nghiệm vẫn là OLS nhưng ổn định số học hơn. Khi X^T X singular, SVD cho pseudo-inverse và coefficient cụ thể không identifiable (nhiều nghiệm cùng minimize loss).

Với 1 feature, model là một đường thẳng trong mặt phẳng \( (x, y) \).

Với 2 feature, model là một mặt phẳng (plane) trong không gian 3D \( (x_1, x_2, y) \). Tưởng tượng tấm bìa nghiêng đặt trong không gian — mỗi điểm dữ liệu \( (x_1, x_2, y) \) được model "ghim" về độ cao trên tấm bìa tại toạ độ \( (x_1, x_2) \).

Với \( n \) feature, model là hyperplane trong không gian \( (n+1) \) chiều — không vẽ được nhưng cùng bản chất: một bề mặt tuyến tính chiều \( n \) chạy xuyên qua đám điểm dữ liệu, sao cho tổng bình phương khoảng cách dọc trục \( y \) giữa điểm và bề mặt là nhỏ nhất.

Trực giác này quan trọng: Linear Regression luôn fit một bề mặt phẳng — không cong. Nếu quan hệ thật là cong (parabola, exponential), bề mặt phẳng sẽ underfit. Đây là motivation cho Polynomial Regression ở Bài 20.

Định nghĩa chuẩn: \( w_j \) = thay đổi kỳ vọng của \( \hat{y} \) khi \( x_j \) tăng 1 unit, các feature khác giữ nguyên (ceteris paribus).

Ví dụ: model giá nhà có \( w_{\text{area}} = 25 \) (triệu đồng / m²), nghĩa là tăng 1 m² thì giá dự đoán tăng 25 triệu — với điều kiện số phòng, vị trí, tuổi nhà không đổi.

Đơn vị của \( w_j \): phụ thuộc unit của \( x_j \) và \( y \). Đổi area từ m² sang ft² → \( w_{\text{area}} \) đổi theo. Đây là lý do không thể so sánh trực tiếp magnitude coefficient giữa các feature có đơn vị khác nhau.

Sau StandardScaler (Bài 8): mọi feature có mean 0, std 1 — đơn vị "1 std". Lúc này coefficient cùng đơn vị, so sánh được:

\[ |w_j^{\text{std}}| \text{ lớn} \Rightarrow x_j \text{ ảnh hưởng } \hat{y} \text{ mạnh hơn (trên thang std)} \]

Cách này được gọi là standardized coefficient hoặc beta coefficient — thường được dùng làm proxy cho feature importance trong Linear Regression.

CẢNH BÁO — coefficient ≠ causation:

• \( w_j \) chỉ phản ánh tương quan có điều kiện (conditional correlation) giữa \( x_j \) và \( y \) trong tập feature đã chọn.
• Thêm/bớt feature khác có thể làm \( w_j \) đổi dấu — gọi là Simpson's paradox.
• Để claim nhân quả, cần thiết kế nghiên cứu phù hợp (RCT, instrumental variable, DAG) — vượt scope ML chuẩn.
• Sai lầm phổ biến: "model fit ra \( w_{\text{cafe}} = 5 \) cho điểm thi → uống cafe giúp tăng điểm". Có thể chỉ vì người chăm học hay uống cafe.

Multicollinearity = hai hoặc nhiều feature có tương quan tuyến tính cao với nhau. Vấn đề chỉ xuất hiện khi có nhiều feature — bài 1 feature không bao giờ gặp.

Hệ quả:

• \( X^T X \) gần singular → nghịch đảo không ổn định.
• Coefficient có variance lớn — thêm/bớt vài sample có thể đảo dấu coefficient.
• Coefficient không identifiable: model có thể chia "đóng góp" giữa 2 feature tương quan theo nhiều cách khác nhau cùng đạt loss tối thiểu.
• Prediction vẫn ổn nếu test set giống train, nhưng diễn giải coefficient không tin được.

Cách phát hiện:

• Correlation matrix — pairwise correlation giữa các feature. Cặp nào \( |r| > 0.8 \) đáng nghi.
• Variance Inflation Factor (VIF) — chuẩn hơn: với mỗi feature \( x_j \), regress nó theo các feature còn lại, lấy \( R^2_j \), rồi: \[ \text{VIF}_j = \frac{1}{1 - R^2_j} \] Quy tắc kinh nghiệm: VIF > 5 đáng lo, VIF > 10 nghiêm trọng. VIF = 1 nghĩa là \( x_j \) hoàn toàn độc lập tuyến tính với các feature khác.

Giải pháp:

• Drop 1 trong 2 feature tương quan cao — chọn cái dễ thu thập / diễn giải hơn.
• Combine feature: PCA (Bài 36) gom các feature tương quan thành thành phần độc lập.
• Regularization: Ridge (Bài 21) ổn định coefficient ngay cả khi multicollinearity; Lasso đẩy bớt coefficient về 0.

Linear Regression chỉ nhận đầu vào số. Feature categorical phải encode trước (Bài 10–11):

• Nominal (không có thứ tự) — One-Hot Encoding.
• Ordinal (có thứ tự) — Ordinal Encoding.

Dummy variable trap: khi OHE với \( k \) category mà giữ đủ \( k \) cột, tổng các cột OHE luôn bằng 1 — trùng với cột bias ngầm. X^T X singular hoàn hảo → OLS lỗi.

Giải pháp: OneHotEncoder(drop="first") hoặc drop="if_binary" để bỏ 1 cột reference. Coefficient của các cột còn lại đọc theo so sánh: "category này so với reference đổi \( \hat{y} \) bao nhiêu".

Lưu ý: với model regularized (Ridge/Lasso) hoặc tree-based, KHÔNG cần drop — vì:

• Tree không quan tâm đến rank của \( X \).
• Ridge thêm \( \lambda I \) vào \( X^T X \) → luôn khả nghịch.

Quy tắc gọn: dùng LinearRegression hoặc OLS thuần → drop cột; các trường hợp khác → tuỳ chọn.

Linear Regression có một điểm dễ gây nhầm: OLS không cần scaling để hoạt động — closed-form cho cùng prediction với mọi phép scale tuyến tính của feature. Coefficient chỉ thay đổi theo tỉ lệ ngược của scale.

Vậy khi nào cần scale?

• So sánh magnitude coefficient — như mục 6, chỉ có ý nghĩa sau khi standardize.
• Khi dùng Gradient Descent (SGDRegressor) — feature chưa scale khiến gradient lệch hướng và hội tụ rất chậm.
• Khi dùng Ridge/Lasso (Bài 21) — penalty \( \lambda \|\mathbf{w}\|^2 \) "phạt" coefficient theo magnitude, nên feature thang lớn (vd thu nhập tính bằng đồng) bị phạt nặng vô lý so với feature thang nhỏ (vd tuổi).
• Khi tính VIF hoặc các diagnostic khác cần ma trận có hệ số ổn định.

Quy tắc thực hành: luôn standardize trừ khi có lý do cụ thể giữ scale gốc. Chi phí scaling rất nhỏ và phòng được nhiều pitfall.

1. Split train/test (Bài 6).
2. Preprocess trong Pipeline (Bài 14):
   • Numeric → StandardScaler.
   • Categorical → OneHotEncoder(drop="first") cho Linear Regression.
3. Fit LinearRegression trên train set.
4. Inspect coefficient sau standardize — sort theo magnitude, kiểm tra dấu có hợp logic không.
5. Diagnostic multicollinearity — correlation matrix + VIF. Nếu có VIF > 10, cân nhắc drop hoặc dùng Ridge.
6. Evaluate trên test set: R² (Bài 18) + RMSE (Bài 17).
7. Inference (optional) — chạy lại bằng statsmodels để lấy p-value, confidence interval (mục 12).

Lưu ý thứ tự: standardize TRƯỚC khi tính VIF, vì VIF đo collinearity theo cấu trúc tương quan, không theo thang đo — kết quả thực tế trùng nhau nhưng quy ước này tránh nhầm lẫn về đơn vị coefficient sau đó.

Dataset California housing có 8 feature số (MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude) và target là giá nhà trung vị theo block (đơn vị 100k USD).

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error

data = fetch_california_housing(as_frame=True)
X, y = data.data, data.target
features = X.columns.tolist()

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("model", LinearRegression()),
])
pipe.fit(X_train, y_train)

y_pred = pipe.predict(X_test)
print(f"R²:   {r2_score(y_test, y_pred):.4f}")
print(f"RMSE: {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")

# In coefficient (đã standardize) sort theo magnitude
coefs = pd.Series(pipe["model"].coef_, index=features)
print("\nCoefficient sort theo |w|:")
print(coefs.reindex(coefs.abs().sort_values(ascending=False).index))

Output mẫu (kết quả có thể khác ±1e-3 tuỳ version sklearn):

R²:   0.5758
RMSE: 0.7456

Coefficient sort theo |w|:
Latitude     -0.8966
Longitude    -0.8689
MedInc        0.8544
AveBedrms     0.3393
AveRooms     -0.3041
HouseAge      0.1226
AveOccup     -0.0396
Population   -0.0023

Đọc: sau standardize, Latitude, Longitude, MedInc là 3 feature có ảnh hưởng mạnh nhất tới giá. Dấu âm của Latitude phù hợp với thực tế California: càng lên Bắc giá thường thấp hơn so với phía Nam (LA, San Diego). Population gần 0 — gần như không đóng góp.

Compute VIF với statsmodels:

from statsmodels.stats.outliers_influence import variance_inflation_factor

X_scaled = StandardScaler().fit_transform(X_train)
vif = pd.DataFrame({
    "feature": features,
    "VIF": [variance_inflation_factor(X_scaled, i) for i in range(X_scaled.shape[1])],
})
print(vif.sort_values("VIF", ascending=False))

Output mẫu:

     feature       VIF
0     MedInc   2.50
2   AveRooms   8.34
3  AveBedrms   6.99
6   Latitude   9.30
7  Longitude   8.94
...

AveRooms và AveBedrms có VIF cao và logic tương quan với nhau (nhà nhiều phòng thường nhiều phòng ngủ). Latitude/Longitude cũng tương quan vì California có hình dạng địa lý đặc biệt. Đây là dấu hiệu cân nhắc Ridge (Bài 21).

Sklearn tập trung vào prediction: fit, predict, score. Không trả về p-value hay confidence interval cho coefficient.

statsmodels bù chỗ thiếu này. Cùng OLS toán học, nhưng output đầy đủ thông tin inference:

import statsmodels.api as sm

# QUAN TRỌNG: phải add_constant để có intercept
X_train_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_train_const).fit()
print(model.summary())

Output cho một bảng dạng:

                            OLS Regression Results
==============================================================================
Dep. Variable:                MedHouseVal   R-squared:                       0.612
Model:                            OLS       Adj. R-squared:                  0.612
...
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -36.94      0.722    -51.17      0.000     -38.36     -35.53
MedInc          0.43      0.005     90.93      0.000       0.42       0.44
HouseAge        0.01      0.001     14.86      0.000       0.01       0.01
...
==============================================================================

Đọc:

• coef — giá trị coefficient (chưa standardize, đơn vị gốc).
• std err — sai số chuẩn của ước lượng coefficient.
• P>|t| — p-value của giả thuyết \( H_0: w_j = 0 \). p < 0.05 → reject H_0, feature "có ý nghĩa thống kê".
• [0.025, 0.975] — confidence interval 95% cho coefficient.

Quy tắc chọn:

• Predict là mục tiêu chính → sklearn.
• Inference (feature có ý nghĩa không, ảnh hưởng bao nhiêu ± khoảng) → statsmodels.
• Cả hai cho cùng coefficient OLS (chỉ khác ở chỗ statsmodels không tự thêm intercept).

• Curse of dimensionality — nhiều feature + ít sample (\( d \) gần \( n \)) → \( X^T X \) gần singular và model overfit. Quy tắc kinh nghiệm: \( n \geq 10d \), lý tưởng \( n \geq 50d \). Nếu không đủ sample, dùng Ridge/Lasso hoặc giảm chiều.
• Multicollinearity bị bỏ qua — predict vẫn cho kết quả "đẹp" trên train, nhưng coefficient không tin được và model fragile trên data mới. Luôn check VIF khi dataset có nhiều feature numeric tương tự.
• Categorical chưa encode — gọi LinearRegression.fit() trên DataFrame có cột string sẽ raise ValueError. Encode trước hoặc bọc trong ColumnTransformer.
• OneHotEncoder không drop cột — gây dummy trap, sklearn không báo lỗi (vì dùng SVD), nhưng coefficient không identifiable.
• Quên add_constant trong statsmodels — model fit không có intercept (qua gốc), R² và coefficient sai. Sklearn tự thêm intercept; statsmodels thì không.
• Diễn giải coefficient sang causation — \( w_j \) chỉ là correlation conditional. Phát biểu kiểu "tăng \( x_j \) sẽ tăng \( y \)" sai về mặt khoa học khi chưa có thiết kế nhân quả.
• Extrapolation — predict với feature value ngoài range train có thể cho kết quả vô lý. Hyperplane kéo dài vô tận không phản ánh ràng buộc thực tế.

Bài 1 — Top 3 feature theo |coef|. Load California housing, split 80/20, fit Pipeline(StandardScaler, LinearRegression). In top 3 feature có \( |w_j| \) lớn nhất. So sánh với output mẫu ở mục 11.

Bài 2 — Correlation matrix. Tính X_train.corr(), hiển thị bằng seaborn.heatmap. Tìm các cặp feature có \( |r| > 0.5 \). Cặp nào trông tự nhiên (vd AveRooms với AveBedrms)?

Bài 3 — Drop feature tương quan cao. Drop một trong cặp AveRooms/AveBedrms, refit pipeline, so sánh R² test với baseline. Drop feature có làm R² giảm nhiều không? Coefficient của các feature còn lại thay đổi thế nào?

Bài 4 — VIF check. Compute VIF cho từng feature (code mục 11). Drop feature có VIF cao nhất, recompute VIF cho tập còn lại — VIF giảm không? Lặp lại đến khi mọi feature có VIF < 5.

Bài 5 — statsmodels inference. Fit lại bằng statsmodels.OLS (nhớ add_constant!). Feature nào có p-value < 0.05? Feature nào CI chứa 0 (nghĩa là coefficient không khác 0 đáng kể về mặt thống kê)?

Gợi ý đáp án Bài 1: Top 3 thường là Latitude, Longitude, MedInc. Dấu của Latitude âm, MedInc dương — phù hợp logic (vùng phía Bắc rẻ hơn, thu nhập cao hơn thì nhà đắt hơn).

Bài 20: Polynomial Regression — model quan hệ phi tuyến — vẫn là Linear Regression, nhưng feature được biến đổi trước qua các bậc \( x, x^2, x^3, \ldots \) (hoặc tương tác \( x_1 x_2 \)). Cho phép fit đường cong mà vẫn dùng OLS chuẩn. Kèm cảnh báo về overfit khi bậc quá cao.