變數選擇

換了問題，挑變數的規則就整套換掉

前面 B2、B6 兩個家族的迴歸都在回答病因學的問題：這個暴露對這個結果有沒有效應。 那種模型裡「要不要放某個變項」是一個因果結構的問題——它是干擾因子（confounder）就必須放， 它是中介（mediator）就不能放，它是對撞因子（collider）放了會製造偏誤。判準在 世代研究與 DAG 那一頁。

預測模型問的是這個人的風險是多少。於是同一個決定換了判準：

三種策略，同一份資料

用 survival::rotterdam（開發，2982 位病人、1713 個事件） 配 survival::gbsg（外部驗證，686 位）。候選變項有 8 個： 年齡、停經狀態、腫瘤大小、分級、淋巴結數、PgR、ER、荷爾蒙治療。

三種在論文裡最常見的做法：

1. 全模型法（full model approach）——事前把 8 個全部指定進去，不論 p 值。
2. 逐步迴歸——用 AIC 從全模型往回刪。
3. 單變量篩選——先一個一個變項各跑一次，p < 0.05 的才進多變量模型。

逐步迴歸的問題不是它挑得不準，是它挑得不穩

模型診斷那一頁已經用 300 次模擬示範過一件事： 把 10 個純亂數變項丟進 stepAIC()，它平均會留下 1.81 個，其中 48.0% 的執行至少出現一個 「顯著」的亂數。那個示範不在這裡重做。

這一頁問一個不同的問題：同一個演算法、同一群人，只是重抽一次樣本，它會給你同一個模型嗎？

把開發世代重抽 200 次，每次都跑一遍相同的向後 AIC，記錄每個變項被留下的比例：

在 200 次重抽樣裡，向後 AIC 一共產生了 25 個不同的模型， 出現最多次的那一個也只佔 22%。 把每次抽樣縮到兩百人（更接近多數單中心研究的規模），變成 69 個不同的模型，最常見的那個只佔 12%。

單變量篩選有一個逐步迴歸沒有的毛病

單變量篩選（先各跑一次、p < 0.05 才進多變量）比逐步迴歸更常見，因為它看起來更保守。它不是。

看上表的停經狀態：單變量 p 值遠小於 0.001，所以它通過篩選、進了多變量模型； 但向後 AIC 把它刪掉了。原因是它和年齡幾乎重疊——校正年齡之後它沒有剩下什麼。

反過來的情況更危險：一個變項可以在單變量下毫無訊號，卻在校正之後變得重要。 單變量篩選會在它有機會被校正之前就把它丟掉，而且丟掉的痕跡不會出現在最終模型的任何一個數字裡。

那到底該怎麼挑

現在的建議（TRIPOD+AI、Riley 等人的方法學文獻）基本上是同一句話： 在看到結果之前就把候選變項與模型形式定下來，然後全部放進去。

動手跑一次

library(survival); library(MASS)
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")
f <- as.formula(paste("Surv(rfs_time, rfs_event) ~", paste(cand, collapse = " + ")))

full <- coxph(f, data = rot)          # 全模型法：就這一行
step <- stepAIC(full, direction = "backward", trace = 0)

# 單變量篩選：注意它問的是「不校正時有沒有關聯」
sapply(cand, function(v)
  summary(coxph(as.formula(paste("Surv(rfs_time, rfs_event) ~", v)),
                data = rot))$coefficients[1, "Pr(>|z|)"])

# 穩定度：同一份資料重抽，看它每次挑到什麼
tally <- setNames(numeric(length(cand)), cand)
for (b in 1:200) {
  d  <- rot[sample(nrow(rot), replace = TRUE), ]
  ff <- f; environment(ff) <- environment()   # 少這行會整批靜默失敗
  st <- stepAIC(coxph(ff, data = d), direction = "backward", trace = 0)
  tally[names(coef(st))] <- tally[names(coef(st))] + 1
}
round(100 * tally / 200)

import numpy as np
import pandas as pd
from lifelines import CoxPHFitter

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"]

# 全模型法
full = CoxPHFitter().fit(d[cols], "rfs_time", "rfs_event")
print(full.concordance_index_)

# 懲罰版本（B5-02 的主題），不是變數選擇的替代品而是它的替代方案
pen = CoxPHFitter(penalizer=0.1, l1_ratio=1.0).fit(d[cols], "rfs_time", "rfs_event")
print(pen.params_[pen.params_ != 0].index.tolist())

當論文用的是機器學習模型時

「XGBoost 預測 AKI」「隨機森林預測敗血症」這類論文在臨床資訊與重症領域已經很多， 而醫學生第一次撞到它們通常是在書報討論上。這一節只教怎麼讀，不教怎麼跑—— 把模型換成森林或梯度提升是另一條技能樹，但驗收標準一格都沒有變。

以下依序談 variable importance、impurity 偏誤、SHAP，以及那件沒有改變的事。

一、variable importance 不是迴歸係數。 它沒有方向、沒有信賴區間、也不是效應量。 一個變項排在第一名，只代表在那個套件所用的重要度定義下它排第一—— permutation importance 的意思是「把它打亂之後模型的預測變差最多」， 而樹系套件的預設往往是下一點要講的 impurity-based 版本，兩者可以排出不同的名次。 不論是哪一種，它都不告訴你風險是升是降、升多少， 更不能被引用成「某因子使風險上升多少」。論文的 Discussion 把重要度排序寫成危險因子清單時， 那是把兩種完全不同的量當成同一個。

二、impurity-based importance 對高基數變項有系統性偏誤。 樹在每個節點挑切點， 取值愈多的變項可切的位置愈多，愈容易在隨機的雜訊上找到一個看起來有用的切法。 所以連續變項與類別很多的變項（例如就診科別、藥品代碼）會系統性地排在前面， 而二元變項會系統性地被低估——這與它們實際上有沒有用無關。 用 permutation importance 或在獨立資料上重算會緩解，但不會消失。

三、SHAP 是局部歸因，不是因果效應。 它回答的是「對這一位病人， 模型的預測值相對於基準被這個變項推了多少」，而那是在模型已經長成這樣的前提下的分解。 模型學到的關聯若來自干擾、選擇或反向因果，SHAP 會忠實地把那份關聯歸因出來， 而且畫得非常有說服力。一張漂亮的 SHAP 圖不是因果證據，它是模型的自述。

四、驗收標準完全沒變。 鑑別力、校準、 決策曲線、 外部驗證——這四項一件都不能少。 而 ML 論文最常跳過的正是校準：樹系模型的預測值天然往中間集中， 校準斜率大於 1 是常見結果，而只報 AUC 的論文不會讓你看到這件事。

讀論文時要問的四件事

1. 候選變項是事前定的，還是資料挑的？ Methods 沒寫「事先指定」就要當作是資料挑的。
2. 挑完之後有沒有重新驗證？ 選變數本身用掉了自由度；事後的 p 值、CI、C-index 全部偏樂觀。
3. 模型裡的變項在部署時真的拿得到嗎？ 這是 data leakage 唯一能被讀者抓到的地方。
4. 有沒有報告被刪掉的那些？ 只列最終模型的論文，讀者無從判斷刪的過程有多勉強。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B5-01-variable-selection.R

延伸觀看

▶What Are Clinical Prediction Rules?

▶AI 臨床研究實戰 EP7｜預測 vs 分類、Data Leakage、Propensity Score

▶Key Steps and Common Pitfalls in Clinical Prediction Model Research

回報這一頁的問題