專家已經雙重審閱，尚未人工抽查

收縮與懲罰迴歸

樣本不夠時，最大概似估出來的係數會太極端——模型把這一批人身上的雜訊當成訊號學了進去。這一頁把同一個模型在七種開發樣本量下各配兩百次，量出過度配適有多大，並示範 Van Houwelingen 的一致收縮因子幾乎可以事先算出「校準斜率會掉到多少」。也講 ridge、lasso、elastic net 的差別，以及為什麼懲罰過的係數不能再當效應量讀。

過度配適不是「模型太複雜」，是「模型把雜訊也記住了」

最大概似估計做的事是：找一組係數，讓手上這批人的資料看起來機率最大。 資料裡有訊號也有雜訊，而這個目標函數不會區分兩者——只要把雜訊也擬合進去能讓概似再高一點，它就會做。

後果很具體：係數的絕對值會偏大。模型把高風險的人推得比實際更高、低風險的人壓得比實際更低。於是在新的人身上，預測值會太極端——這在校準那一頁會表現為 校準斜率小於 1。

這件事可以直接做出來。survival::rotterdam 有 2982 位病人，夠大到可以製造一個過度配適的情境：每次只抽一小部分人來配模型，再把模型丟到剩下的人身上量表現。剩下的人來自同一個世代、同一批測量，所以兩者的差距只可能是過度配適，不會混到族群差異。

同一個 8 變項的 Cox 模型，七種開發樣本量，每種各重複 200 次：

兩張折線圖。左圖橫軸是開發樣本量（對數尺度），兩條線分別是在同一批資料上量到的表面 C-index 與在剩下病人身上量到的 C-index；樣本量小的時候兩條線相差很遠，隨樣本量增加逐漸收斂。右圖是啟發式收縮因子與保留樣本上的校準斜率，兩條線幾乎重疊，都從約 0.55 上升到接近 1。 — 左：表面表現（在配模型的同一批資料上量）永遠比保留樣本上的表現好，而且樣本愈小差距愈大。右：Van Houwelingen 的啟發式收縮因子（只用開發資料就能算）與實際在保留樣本上量到的校準斜率，兩條線幾乎重合。產圖腳本 `figures/scripts/B5-02-shrinkage.R`

開發樣本量	事件數	EPV	表面 C-index	保留樣本 C-index	樂觀偏誤	保留樣本校準斜率
80	46	5.7	0.697	0.628	0.0688	0.527
120	69	8.6	0.683	0.639	0.0438	0.648
200	115	14.4	0.679	0.651	0.0283	0.779
350	201	25.1	0.674	0.657	0.0177	0.852
600	345	43.1	0.672	0.662	0.0102	0.909
1200	689	86.2	0.670	0.665	0.0048	0.960
2400	1379	172.3	0.670	0.668	0.0019	0.994

收縮：把係數往零拉一點

既然問題是「係數太極端」，最直接的修法就是把它們往零拉。這就是收縮（shrinkage）。

拉多少？Van Houwelingen 與 Le Cessie 給了一個只用開發資料就能算的啟發式數字：

\hat{s} = \frac{\chi^2_{\mathrm{LR}} - p}{\chi^2_{\mathrm{LR}}}

分子是模型的概似比卡方值減掉參數個數，分母是卡方值本身。把每一個係數都乘上這個 一致收縮因子（uniform shrinkage factor），模型就不再那麼極端。

它為什麼合理：卡方值裡有一部分是真訊號、有一部分是「多估 p 個參數本來就會多出來的」，後者的期望值恰好是 p。把後者扣掉再取比例，得到的就是「訊號佔多少」。

這個數字能不能事先算出實際會掉到多少？ 上面那張圖的右半邊就是在檢查這件事—— 啟發式收縮因子（只看開發資料）與實際在保留樣本上量到的校準斜率，在七種樣本量下都貼得很近。在最小的樣本量下它預測 0.591，實際量到 0.527。

ridge、lasso、elastic net：同一個想法的三種罰則

懲罰迴歸把收縮直接寫進估計過程：不再單純最大化對數概似，而是最大化

\ell(\beta) - \lambda \left[ \frac{1-\alpha}{2}\sum_j \beta_j^2 + \alpha \sum_j |\beta_j| \right]

括號裡就是罰則， $\lambda$ 決定罰多重、 $\alpha$ 決定罰的形狀：

	罰則	會不會把係數變成剛好 0	典型用途
ridge	係數平方和（alpha 為 0）	不會，只會全部縮小	所有候選變項都想留、只是想讓它們別太極端
lasso	係數絕對值和（alpha 為 1）	會，等於同時做變數選擇	候選變項多、想要一個比較短的模型
elastic net	兩者的混合	會，但比 lasso 溫和	變項之間高度相關時（lasso 會在一群相關變項裡隨機挑一個留）

差別的來源是罰則在零點的形狀：絕對值在零點有一個尖角，最佳解會被「卡」在零上；平方在零點是平滑的，所以只會逼近零而不會到零。

拿 rotterdam 抽出 200 位病人（125 個事件，EPV 15.6）當開發資料，其餘 2782 位當保留樣本：

左圖是 lasso 的係數路徑，橫軸是懲罰強度的對數，八條線代表八個預測變項乘上各自的標準差後的係數，隨著懲罰增強逐一被壓到剛好零；一條垂直虛線標出交叉驗證選出的懲罰值。右圖是同一份資料重跑六十次交叉驗證所選出的懲罰值直方圖，分布明顯右偏，最大值約為最小值的三點七倍。 — 左：lasso 路徑。係數乘上各自預測變項的標準差之後才畫在同一張圖上，否則以年、毫米、0/1 為單位的變項不能比較。垂直虛線是交叉驗證選出的懲罰強度。右：同一份資料、同一段程式，只換交叉驗證的隨機分組，重跑六十次得到的懲罰值。產圖腳本 `figures/scripts/B5-02-shrinkage.R`

做法	非零係數個數	保留樣本 C-index	保留樣本校準斜率
最大概似（完全不加懲罰）	8	0.6419	0.807
一致收縮因子	8	0.6419	0.962
ridge（alpha = 0）	8	0.6532	1.268
elastic net	5	0.6517	1.010
lasso（alpha = 1）	5	0.6501	0.953

懲罰過的係數不是效應量

lasso 在這份資料上把 3 個變項壓成零。看一下四種做法的係數：

變項	最大概似	一致收縮後	ridge	lasso
age	0.019	0.016	0.005	0.007
menopausal	-0.183	-0.153	0.022	0
tumour size	0.009	0.008	0.009	0.009
grade 3	0.639	0.536	0.359	0.517
nodes	0.071	0.059	0.048	0.067
log PgR	0.006	0.005	-0.015	0
log ER	-0.115	-0.097	-0.054	-0.087
hormonal tx	-0.062	-0.052	0.008	0

這張表的每一個非零數字都不能寫成「校正後的 log hazard ratio」。

理由有三層：

懲罰估計是刻意有偏的。它用偏誤換變異，整體預測會比較好，但單一係數的期望值不再是真值。
它沒有可用的標準誤。glmnet 不會給你 p 值與信賴區間，這不是功能缺漏而是刻意的—— 懲罰之後的抽樣分布沒有簡單的形式，硬套常態近似會給出過窄的區間。
lasso 的「零」不代表「無關」。它代表在這個 $\lambda$ 、這一批人身上，那個變項的貢獻小到不值得付罰則的代價。換一批人、換一個 $\lambda$ ，它可能就回來了。

所以懲罰模型的產出只有一種正當用法：把整組係數當成一台計算風險的機器。要談某個變項的效應，回去配一個不加懲罰的模型，並且準備好回答干擾因子那一套問題。

懲罰參數本身也是估計出來的

$\lambda$ 通常用交叉驗證選：切成 k 份，對一整排候選 $\lambda$ 各算一次交叉驗證的偏差，挑最小的那個。

問題是分組本身是隨機的。同一份資料、同一段程式，只換隨機分組重跑 60 次：

	最小	中位數	最大
交叉驗證選出的 lambda	0.0166	0.0264	0.0610
留下的非零係數個數	5	5	5

最大值是最小值的 3.68 倍。

收縮不能救的東西

懲罰不是樣本量不足的解藥。 上面那張表已經給了提示：收縮把校準斜率拉回來了，但這件事本身不會憑空生出鑑別力。分開看兩件事：一致收縮把每個係數乘上同一個因子，不改變任何人的相對排序，保留樣本的 C-index 與最大概似完全相同 （都是 0.6419）——這是數學上的必然，不是量出來的差距。 ridge 對不同係數的收縮程度不同，所以確實會動到排序：這一次的點估計從 0.6419 到 0.6532，高出 0.0113，是這張表裡最大的一組差距；但那是單次切分的點估計，本頁沒有估計它的不確定性。

收縮做的是「承認我知道得沒那麼多」，它讓模型的預測值誠實一點，不會讓資料裡本來沒有的訊號出現。當樣本量嚴重不足時，方法學文獻的建議是乾脆不要開發這個模型，而不是靠懲罰硬撐。判斷樣本夠不夠的工具在下一頁。

動手跑一次

library(survival); library(glmnet)
data(cancer, package = "survival")

rot <- rotterdam
rot$rfs_time  <- pmin(rot$rtime, rot$dtime)
rot$rfs_event <- as.integer(rot$recur == 1 | rot$death == 1)
rot$size_mm   <- c("<=20" = 15, "20-50" = 35, ">50" = 60)[as.character(rot$size)]
rot$grade3    <- as.integer(rot$grade >= 3)
rot$log_pgr   <- log1p(rot$pgr); rot$log_er <- log1p(rot$er)
cand <- c("age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon")

set.seed(1)
i     <- sample(nrow(rot), 200)          # 刻意做出一個樣本不足的情境
small <- rot[i, ]; rest <- rot[-i, ]

fit <- coxph(as.formula(paste("Surv(rfs_time, rfs_event) ~",
                              paste(cand, collapse = " + "))), data = small)

# 一致收縮因子：只要配好的模型就能算
lr <- 2 * diff(fit$loglik); p <- length(coef(fit))
(shrink <- (lr - p) / lr)
coef(fit) * shrink                       # 收縮後的係數

# 懲罰迴歸
x <- as.matrix(small[, cand]); y <- Surv(small$rfs_time, small$rfs_event)
cv_lasso <- cv.glmnet(x, y, family = "cox", alpha = 1, cox.ties = "breslow")
coef(cv_lasso, s = "lambda.min")         # 被壓成 . 的就是被挑掉的

# 懲罰參數的穩定度：重跑就知道
lam <- replicate(60, cv.glmnet(x, y, family = "cox", alpha = 1,
                               cox.ties = "breslow")$lambda.min)
range(lam); max(lam) / min(lam)

# 在保留樣本上驗收：C-index 與校準斜率要一起看
b  <- as.numeric(coef(cv_lasso, s = "lambda.min"))
lp <- as.matrix(rest[, cand]) %*% b
concordance(Surv(rest$rfs_time, rest$rfs_event) ~ lp, reverse = TRUE)$concordance
coef(coxph(Surv(rfs_time, rfs_event) ~ lp, data = rest))

驗證環境：R 4.6.0 + survival 3.8.6 + glmnet 5.0。glmnet 的 Cox 模型預設用 Breslow 處理同時間事件，與 coxph 的 Efron 不同；把 cox.ties 明確寫出來可以避免跨版本漂移。glmnet 內部會標準化預測變項，但回傳的係數已經換算回原始尺度。

import numpy as np
import pandas as pd
from lifelines import CoxPHFitter
from lifelines.utils import concordance_index

RD = "https://vincentarelbundock.github.io/Rdatasets/csv/"
d = pd.read_csv(RD + "survival/rotterdam.csv")

# 與上面 R 端相同的衍生欄位。rotterdam 原始欄位沒有這些，少了這一段
# 下面每一行都會 KeyError。
d["rfs_time"]  = d[["rtime", "dtime"]].min(axis=1)
d["rfs_event"] = ((d["recur"] == 1) | (d["death"] == 1)).astype(int)
d["size_mm"]   = d["size"].map({"<=20": 15, "20-50": 35, ">50": 60})
d["grade3"]    = (d["grade"] >= 3).astype(int)
d["log_pgr"]   = np.log1p(d["pgr"]); d["log_er"] = np.log1p(d["er"])
cand = ["age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon"]
cols = cand + ["rfs_time", "rfs_event"]

rng = np.random.default_rng(1)
i = rng.choice(len(d), 200, replace=False)
small, rest = d.iloc[i], d.drop(d.index[i])

ml    = CoxPHFitter().fit(small[cols], "rfs_time", "rfs_event")
lasso = CoxPHFitter(penalizer=0.05, l1_ratio=1.0).fit(small[cols], "rfs_time", "rfs_event")
ridge = CoxPHFitter(penalizer=0.05, l1_ratio=0.0).fit(small[cols], "rfs_time", "rfs_event")

for name, m in [("ml", ml), ("lasso", lasso), ("ridge", ridge)]:
    lp = rest[cand].to_numpy() @ m.params_[cand].to_numpy()
    # 注意 concordance_index 的方向：風險分數要取負號才是「存活時間」的順序
    print(name, concordance_index(rest["rfs_time"], -lp, rest["rfs_event"]))

scikit-learn 沒有懲罰 Cox；lifelines 的 CoxPHFitter 用 penalizer 與 l1_ratio 兩個參數對應 glmnet 的 lambda 與 alpha，但兩者的 lambda 尺度定義不同，不要直接搬數字。

常見誤用

誤用	為什麼錯
把 lasso 選出來的變項當成「重要的預測因子」	那是在某個懲罰強度下的取捨結果，換一批人就會變
報告懲罰後係數的 p 值或信賴區間	懲罰估計沒有簡單的抽樣分布，常態近似的區間會太窄
用懲罰迴歸解決樣本量不足	它讓預測值誠實，不會變出資料裡沒有的訊號
期待收縮讓 AUC 上升	一致收縮完全不改排序，C-index 一位小數都不會動
只跑一次 `cv.glmnet` 就定案	分組是隨機的，重跑會給出不同的懲罰強度
把懲罰強度的選擇留在重抽樣迴圈外	那也是從資料裡做的決定，樂觀偏誤照樣進來
認為 ridge 一定比不加懲罰保險	交叉驗證選出的懲罰可能過重，校準斜率會跑到 1 以上
在沒有標準化的情況下比較懲罰模型的係數大小	罰則對所有係數一視同仁，尺度大的變項會被不成比例地罰
用懲罰模型的係數談因果	懲罰是為了預測而刻意引入偏誤的，不是效應估計

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B5-02-shrinkage.R

讀讀看這張圖

答案取自產生本頁圖表的同一份統計輸出，不是另外打上去的。

同一個八變項的 Cox 模型，在只有八十人的開發樣本上表面 C-index 是 0.697，把模型丟到剩下的人身上只剩 0.628。這個落差說明什麼？

看答案與解析

正確答案: 那個落差就是 0.069 的樂觀偏誤，而剩下的人來自同一個世代、同一批測量，所以它不可能是族群差異

0.069 就是表面表現減掉保留樣本表現。這裡刻意讓保留樣本來自同一個世代，正是為了把族群差異這個解釋關掉，差距只剩過度配適一個來源。0.002 是樣本量最大那一列的樂觀偏誤，它恰恰證明了相反的事：同一份資料、同一個模型，只是樣本變大，過度配適就幾乎消失，所以那不是資料的屬性，是模型複雜度與樣本量之間的關係。0.527 那個校準斜率也不是無關的第三件事——係數太極端正是預測值太極端的原因，鑑別力的落差與校準斜率掉下來是同一個機制的兩個面向。

一致收縮因子把每一個係數乘上同一個常數。表上它的保留樣本 C-index 與最大概似完全一樣，校準斜率卻從 0.807 拉到 0.962。為什麼？

看答案與解析

正確答案: C-index 只看排序，而所有係數乘上同一個正數不改變排序，所以它與最大概似同為 0.6419 是必然的

0.6419 出現兩次不是巧合：一致收縮把每個係數乘上同一個因子，任兩個人的線性預測值大小關係完全不變，而 C-index 只問排序對不對。ridge 的 0.6532 與 lasso 的 0.6501 之所以會動，是因為它們對不同係數的收縮程度不同——那才是排序改變的原因，不是「收縮得比較徹底」。這一題真正要記住的是分工：收縮修的是校準（斜率從 0.807 拉到 0.962），不是鑑別力，期待懲罰讓 AUC 上升多半會失望。

同一份保留樣本上，ridge 的校準斜率是 1.268。這個數字該怎麼讀？

看答案與解析

正確答案: 1.268 大於 1 代表縮過頭了：預測值變得太保守、風險差異被壓得太小，所以懲罰不是愈重愈好

校準斜率的目標是 1，兩個方向的偏離都是問題。0.962 離 1 只有幾個千分點，1.268 離 1 有將近三成，說兩者差不多只要把數字擺在一起就站不住。方向也不同：小於 1 是預測值太極端，大於 1 是預測值太保守、把高低風險壓在一起。所以 ridge 不是「修得比較徹底」，它是修過了頭；最大概似的 0.807 與 ridge 的 1.268 是兩種相反的錯，不是同一條路上的兩個位置。交叉驗證選出的懲罰強度是為了讓某個預測誤差最小，它從來沒有以「校準斜率等於 1」為目標。

lasso 把停經狀態的係數壓成剛好零，而同一個變項在最大概似下的係數是 -0.183。哪一個說法對？

看答案與解析

正確答案: -0.183 是不加懲罰時的估計；lasso 的零只代表在這個懲罰強度、這一批人身上，這個變項不值得付罰則的代價

lasso 的零是一個最佳化的結果，不是一個檢定的結論。0.022 與 -0.183 方向相反確實刺眼，但那不是誰算錯：ridge 在相關變項之間分配權重的方式與最大概似不同，係數的方向在懲罰之後本來就不必保持一致，這正是「懲罰後的係數不能當效應量讀」的其中一個理由。一致收縮的 -0.153 則是把原來的估計乘上同一個因子，它與 lasso 的零不是同一種操作：前者不改變任何變項的相對權重，後者是在罰則的尖角上把某些變項卡在零。真的要談這個變項有沒有效應，回去配一個不加懲罰的模型，並且準備好回答干擾因子那一套問題。

同一份資料、同一段程式，只換交叉驗證的隨機分組重跑六十次，選出來的懲罰強度最小 0.0166、最大 0.0610，而留下的非零係數個數一次都沒變。這該怎麼讀？

看答案與解析

正確答案: 中位數 0.0264 才是比較可信的落點；變項個數沒變是因為這份資料的強弱變項分得很開，中間沒有邊緣變項

0.0264 是六十次重跑的中位數，實務上的建議正是取它，或改用比較穩定的 lambda.1se，而不是相信單次執行給的 0.0219——一個數字落在自己那個分布的範圍裡，本來就不能當成它可信的理由，何況那個範圍從 0.0166 到 0.0610。變項個數沒有變也不是穩定的證據：這份資料裡淋巴結數、分級、腫瘤大小遠遠強過其他變項，整段懲罰區間都還壓不掉它們，而弱的那些早就被壓掉了，中間沒有處在邊緣的變項可以翻面。3.6784 倍的晃動換一份強弱沒分那麼開的資料就會改變模型組成。最後一點更重要：選懲罰強度本身也是從資料裡做出來的決定，它必須一起進到重抽樣迴圈裡。

Van Houwelingen 的一致收縮因子只用開發資料就能算。在最小的那個樣本量下它給出 0.591，而實際在保留樣本上量到的校準斜率是 0.527。這個對照的意義是什麼？

看答案與解析

正確答案: 0.591 是配完模型當場就能算出來的，而它跟事後才量得到的校準斜率貼得很近；它的價值在於便宜

0.591 與 0.527 貼得近，而且在七種樣本量下都貼得近，這正是本頁那張圖右半邊在檢查的事。0.980 那一列不能拿來說「只在大樣本才準」——那一列的校準斜率同樣接近 1，兩者一樣是貼近的，它其實是「樣本充足時過度配適很小」這個事實的另一種寫法。0.873 與同一列的校準斜率不完全相等也不是巧合的證據：一致收縮因子是啟發式的，它假設所有係數都該按同一個比例縮，而真實情況通常不是，所以它給的是一個很好的近似而不是一個等式。要注意這份吻合是這份資料上的觀察，不是保證。

用到這個方法的章節

預測模型研究

延伸觀看

ENEhsan Karim· 30 min中段講過度配適與它對校準的影響，正是本頁第一、二節的內容。

ENNIHR Maudsley BRC· 62 min一小時的建模全景，懲罰迴歸放在整個流程裡的位置在這裡看得比較清楚。

素材來源與授權

本頁為原創內容