Model Comparison / Feature Selection

1. Big model is not necessarily good model（slides 1–3）

呢講用 multiple linear regression 做例子，講點樣喺「大 model」同「細 model」之間揀。Slide 1 列出成講嘅路線：

• Discouraging big model by penalty term（$C_p$、AIC、BIC）
• 呢啲諗法用嚟做 variable selection（feature selection）
• Prediction ability vs cross validation

Module 4 講過 OLS 嘅「best result」係 minimize Residual sum of squares (RSS)：

點樣 compare 大細 model？本講提供兩個方法：

	Method 1：Penalty（slides 3–6）	Method 2：Cross Validation（slides 7–12）
諗法	Punish big model：RSS + penalty term	直接睇 prediction ability on unseen data
工具	Mallows' $C_p$、AIC、BIC	Train-test split、K-fold CV、LOOCV
用嘅 data	Full training set（fit 一次）	Split data，fit 多次
準則	Lower is better	CV error（MSE）lower is better

2. Method 1：用 Penalty 懲罰大 model（slides 3–6）

2.1 Mallows' $C_p$（slide 4）

點理解：$RSS_{sub}$ 隨 $d$ 增加而跌（fit 得好啲），penalty $2d\hat\sigma^2$ 隨 $d$ 增加而升。$C_p$ 係兩股力嘅 trade-off：新加嘅 predictor 要令 RSS 跌多過 $2\hat\sigma^2$，先值得加。

2.2 AIC 同 BIC（slide 5）

除咗 $C_p$，仲有其他方法懲罰大 model：

Criterion	每個 predictor 嘅 penalty	特點
$C_p$	$2\hat\sigma^2$（再除 $n$）	Lower is better
AIC	$2\hat\sigma^2$（再除 $n\hat\sigma^2$）	$AIC=C_p/\hat\sigma^2$，同 $C_p$ 成正比
BIC	$(\log n)\,\hat\sigma^2$	Penalty 最重 → 揀較細 model

2.3 Worked example

$n=50$，用 $\hat\sigma^2=13.33$（全部 models 共用）。比較兩個 models：

Model	$d$（含 intercept）	$RSS_{sub}$
A（2 個 predictors）	3	650
B（4 個 predictors）	5	600

B 嘅 RSS 細啲（意料之內，因為 model 大啲）。但計埋 penalty：

• $C_p(A)=\frac{1}{50}\left[650+2(3)(13.33)\right]=\frac{730}{50}=14.60$
• $C_p(B)=\frac{1}{50}\left[600+2(5)(13.33)\right]=\frac{733.3}{50}=14.67$
• $BIC(A)=\frac{1}{50(13.33)}\left[650+(\log 50)(3)(13.33)\right]=\frac{650+156.5}{666.7}\approx1.210$
• $BIC(B)=\frac{1}{666.7}\left[600+(\log 50)(5)(13.33)\right]=\frac{600+260.8}{666.7}\approx1.291$

兩個 criteria 都係 A 較低 → 揀 A。B 多咗 2 個 predictors，RSS 只跌咗 50，唔夠抵銷 penalty。留意 BIC 下 A 贏得更多，因為 BIC penalty 更重。

Slide 6：呢啲定義背後嚴格嘅數學 rationale 涉及更高深嘅統計知識，屬 Out-of-syllabus。

3. Method 2：Cross Validation（slides 7–12）

Crucial idea 一樣係 big model is not necessarily good model，但 Method 2 直接問：個 model 喺 unseen data 上係咪真係 work？方法係比較大細 models 嘅 prediction abilities。Underlying principle：Train-test split。

3.1 Train-test split（slides 8–9）

（$m$ = test data 數目。）Slide 9 嘅圖：data 分成 Training data 同 Validation/Test data 兩段。

簡單 split 有個弱點：結果視乎邊啲 data 啱啱好分咗落 test set，換一個 split 可能得出唔同 MSE。改良版本喺實際好常用：K-fold CV 同 Leave-one-out。

3.2 K-fold cross validation（slide 10）

Slide 10 嘅圖（5-fold CV，取自 ISLR Figure 5.5）：將 n 個 observations randomly 分成 5 個 non-overlapping groups。每一份輪流做 validation set（test data），其餘 4 份做 training set。Test error 用 5 個 MSE 嘅平均估計。

3.3 Leave-one-out（slides 11–12）

方法	Fit 幾多次	特點
Train-test split	1	最快；結果受 split 方式影響
K-fold CV（常用 K = 5 / 10）	K	平衡準確度同計算量
LOOCV	n	= n-fold CV；最貴

4. Best subset selection（slides 13–17）

（Slide 14 只係一張 savings jar 嘅圖，引入下面嘅 savings dataset。）

4.1 有幾多個 models？（slide 15）

例子：Y = sr（savings rate），4 個 predictors：pop15、pop75、dpi、ddpi。可能嘅 linear models 例如 sr ~ 1、sr ~ pop15+dpi、sr ~ pop75+dpi+ddpi……

4.2 Python：ExhaustiveFeatureSelector + CV（slide 16）

import pandas as pd
import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from mlxtend.feature_selection import ExhaustiveFeatureSelector as EFS

savings = pd.read_csv("savings.csv")
X = savings[['pop15', 'pop75', 'dpi', 'ddpi']]
X = sm.add_constant(X)          # 加一個全 1 嘅 column 代表 constant term（column 0）
y = savings['sr']

efs = EFS(LinearRegression(),
          min_features=1,
          max_features=4,
          fixed_features=(0,),       # column 0（const）一定要包括
          scoring='neg_mean_squared_error',  # higher is better，所以用負 MSE
          cv=5)                      # 5-fold CV
efs = efs.fit(X, y)
print('Best subset:', efs.best_feature_names_)
efs.subsets_                         # 每個 subset 嘅詳細結果
# Something wrong, why cannot cover all 16 possible choices?

Argument	意思
min_features / max_features	Subset 最少 / 最多有幾多個 columns（const 都計一個）
fixed_features=(0,)	第 0 個 column（const）每個 subset 都包括
scoring='neg_mean_squared_error'	mlxtend / sklearn 嘅 score 係 higher is better，所以將 MSE 變負數
cv=5	用 5-fold CV 計每個 subset 嘅 score

4.3 Python：用 BIC 做 subset selection（slide 17）

# 1. 自定 BIC scorer
def bic_scorer(estimator, X, y):
    n_samples = X.shape[0]
    n_features = X.shape[1]
    y_pred = estimator.predict(X)
    rss = np.sum((y - y_pred) ** 2)          # Residual Sum of Squares
    if rss < 1e-10:                          # 避免 log(0)
        return -np.inf
    bic = n_samples * np.log(rss / n_samples) + n_features * np.log(n_samples)
    return -bic                              # mlxtend 係 maximize score，所以回傳負 BIC

# 2. Base estimator
lr = LinearRegression()

# 3. ExhaustiveFeatureSelector
efs = EFS(estimator=lr,
          min_features=1,
          max_features=5,                    # const + 4 個 predictors，16 個 models 齊晒
          fixed_features=(0,),
          scoring=bic_scorer,                # 直接傳 function
          cv=0,                              # 0 = 用成個 training set 評估，唔做 CV
          print_progress=True)
efs = efs.fit(X, y)
print('Best subset:', efs.best_feature_names_)
efs.subsets_

5. Forward and Backward Selection（slides 18–22）

Best subset 對大 p 唔實際（impractical for large p）。Slide 18 問：p = 20 時有幾多個 models？

如果再用 5-fold CV，就要 fit 超過 500 萬次。Forward 同 backward selection 係 sub-optimal but practical 嘅方法。

計算量比較（p = 20）：forward selection 第 1 步試 20 個、第 2 步 19 個……最多 $1+20+19+\cdots+1=1+\frac{20\times21}{2}=211$ 個 models，遠少過 1,048,576。

5.1 Python：SequentialFeatureSelector（slides 20–22）

用 mlxtend.feature_selection.SequentialFeatureSelector，以 forward=True / forward=False 切換。

from mlxtend.feature_selection import SequentialFeatureSelector as SFS

sfs = SFS(LinearRegression(),
          k_features="best",                 # 自動揀 score 最好嘅 subset size
          fixed_features=(0,),
          forward=True,                      # Forward selection；False = Backward
          scoring='neg_mean_squared_error',
          cv=5)
sfs = sfs.fit(X, y)
print('Best MSE score:', -sfs.k_score_)      # 負負得正，變返 MSE
print('Best subset:', sfs.k_feature_names_)
sfs.subsets_

Slide 21（forward）同 slide 22（backward）嘅 code 唯一分別係 forward=True / forward=False。留意呢度用嘅係 CV MSE 做準則，唔係 $C_p$/AIC/BIC。

6. Retrieve the optimal model：Pipeline（slide 23）

from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error

estimator = LinearRegression()
pipeline = Pipeline([('feature_selection', sfs),
                     ('regression', estimator)])

pipeline.fit(X_train, y_train)          # 喺 training data 上做 feature selection + fit
pre_val = pipeline.predict(X_test)      # 直接畀未揀 features 嘅 test data，pipeline 自己揀返 columns

act_val = y_test
bias = np.mean(pre_val - act_val)
rmse = np.sqrt(mean_squared_error(act_val, pre_val))
print(f'Out-of-sample RMSE: {rmse}')

自我測驗 Self-test

Q1. Model A 有 3 個 predictors，Model B 喺 A 嘅基礎上多加 2 個。邊個 training RSS 較細？係咪代表 B 較好？

B 嘅 RSS 一定 ≤ A（bigger model, smaller RSS）。但唔代表 B 較好：small RSS 只係 in-sample fit 好，唔保證 out-of-sample prediction 好，B 可能 overfit。要用 $C_p$/AIC/BIC 或 CV 比較。

Q2. 寫出 Mallows' $C_p$，指出 penalty term，並解釋 $d$ 嘅意思。

$C_p=\frac1n\left[RSS_{sub}+2d\hat\sigma^2\right]$。Penalty term 係 $2d\hat\sigma^2$。$d$ 係 sub-model 嘅 predictors 數目（包括 intercept）。Lower is better。

Q3. $n=100$，$\hat\sigma^2=4$。Model A：$d=3$，RSS = 380；Model B：$d=6$，RSS = 350。用 $C_p$ 同 BIC 揀 model。

$C_p(A)=\frac{380+2(3)(4)}{100}=4.04$；$C_p(B)=\frac{350+2(6)(4)}{100}=3.98$ → $C_p$ 揀 B。
$\log 100\approx4.605$。$BIC(A)=\frac{380+4.605(3)(4)}{400}=\frac{435.3}{400}\approx1.088$；$BIC(B)=\frac{350+4.605(6)(4)}{400}=\frac{460.5}{400}\approx1.151$ → BIC 揀 A。
呢個正好示範 BIC penalty 較重，揀較細 model。

Q4. 幾時 BIC 嘅 penalty 會重過 AIC？

BIC 每個 predictor 嘅系數係 $\log n$，AIC 係 2。$\log n>2\iff n>e^2\approx7.39$，即 $n\ge8$ 時 BIC 已經較重。

Q5. 描述 5-fold CV，並寫出 CV error。

將 data 隨機分 5 份；每次用 1 份做 test、4 份做 training，計 MSE；重複 5 次令每份都做過 test；$CV_{(5)}=\frac15\sum_{k=1}^{5}MSE_k$。

Q6. 有 500 個 data instances，LOOCV 要 fit 幾多次 model？佢同 K-fold CV 有咩關係？

500 次。LOOCV = n-fold CV（K = n），每次留 1 個做 test，所以 computationally expensive。

Q7. 有 6 個 predictors，best subset selection 要考慮幾多個 models？p = 20 呢？

$2^6=64$（包括 null model）。$2^{20}=1{,}048{,}576$，所以大 p 時要用 forward / backward selection。

Q8. Slide 16 嘅 EFS 用 max_features=4，點解唔夠 16 個 models？點改？

X 加咗 const 之後有 5 個 columns，const 係 fixed 而且會計入 feature 數目。max_features=4 即最多 3 個真 predictors，full model 被排除，只得 15 個。改做 max_features=5。

Q9. 簡述 forward selection 三個步驟。佢一定搵到 best subset 嗎？

(1) 由 null model 開始；(2) 試每個餘下 predictor，揀改善 RSS 最多嘅加入；(3) 用 $C_p$/AIC/BIC 比較新舊 model，冇改善就停，否則返 Step 2。唔保證搵到 best subset（greedy、sub-optimal），但計算量細好多。

Q10. 點解 scoring='neg_mean_squared_error' 要用負數？BIC scorer 點解 return -bic？

mlxtend / sklearn 嘅 selector 永遠 maximize score（higher is better）。MSE 同 BIC 都係 lower is better，加負號之後「最大」就等於原本「最細」。

Q11. Slide 23 點解要用 Pipeline？

將 feature selection 同 regression 串埋：只用 training data 揀 features 同 fit model；predict 時直接畀原始 test data，pipeline 自動揀返同一批 columns。咁 out-of-sample RMSE 先係公平嘅評估，避免 test data 資訊漏入 selection。

一頁總結 Cheat sheet