缺失值與多重插補
完整個案分析不是中性的預設，它是一個有假設的選擇——本章用一份天然有缺失的資料證明被排除的人與留下的人確實不同，並把多重插補拆開來寫，讓你看到它其實是「迴歸加雜訊」加上兩行算術。
「完整個案分析」是一個選擇，不是預設
coxph()、lm()、glm() 遇到有缺失的列會安靜地跳過。不警告、不報錯，只是模型的 n 比你以為的小。
這個行為叫完整個案分析（complete-case analysis），它隱含一個很強的假設： 被排除的人與留下的人，在與分析有關的方面沒有系統性差異。
這一頁要做的第一件事，就是檢查那個假設在實際資料上成不成立。
這份資料的缺失長什麼樣
survival::pbc 有 418 人、20 個變項，其中
12 個變項有缺失。缺得最多的幾個：
| 變項 | 缺失數 | 缺失比例 |
|---|---|---|
trig | 136 | 32.5% |
chol | 134 | 32.1% |
copper | 108 | 25.8% |
trt | 106 | 25.4% |
ascites | 106 | 25.4% |
hepato | 106 | 25.4% |
spiders | 106 | 25.4% |
alk.phos | 106 | 25.4% |
figures/scripts/B8-01-missing-data.R那一整塊是重點。 這份資料裡有 312 人進入試驗、106 人拒絕隨機分配 而進入登錄追蹤。試驗檢驗是 protocol 的一部分，所以登錄組依定義沒有那些值：
| 檢驗 | 試驗組缺失 | 登錄組缺失 |
|---|---|---|
chol | 28 / 312（9%） | 106 / 106（100%） |
copper | 2 / 312（1%） | 106 / 106（100%） |
trig | 30 / 312（10%） | 106 / 106（100%） |
platelet | 4 / 312（1%） | 7 / 106（7%） |
alk.phos | 0 / 312（0%） | 106 / 106（100%） |
ast | 0 / 312（0%） | 106 / 106（100%） |
platelet 是這張表上的例外：它兩組都只有零星缺失（登錄組 7 / 106），
不是「登錄組依定義沒有」的結構性缺失，而是一般的零星缺失。chol 與 trig 在試驗組也有非零缺失。
真正 100% 落在登錄組的是 alk.phos 與 ast。
缺失不是隨機灑在資料上的，它有結構，而且那個結構與病人的一個實質特徵（願不願意進入試驗）綁在一起。
三種缺失機制
| 機制 | 意思 | 完整個案分析會怎樣 |
|---|---|---|
| MCAR（完全隨機缺失） | 缺失與任何變項都無關，包括缺失值本身 | 不偏，只是損失效率 |
| MAR（隨機缺失） | 缺失可以由已觀測到的變項解釋 | 不一定有偏：只要在模型共變項給定之後，缺失與結果無關，迴歸係數仍然不偏，只損失效率；這個條件不成立時才有偏，而平均數這類邊際量即使條件成立也可能有偏，因為完整個案是被共變項篩選過的一群人。多重插補兩者都能處理 |
| MNAR（非隨機缺失） | 缺失取決於沒觀測到的值本身 | 有偏；插補救不了，需要敏感度分析 |
被排除的人跟留下的人不一樣
模型用五個變項，完整個案剩 310 人， 排除了 108 人（25.8%）。
那被排除的人長什麼樣？
| 變項 | 留下（n = 310） | 排除（n = 108） | SMD |
|---|---|---|---|
| age | 49.95 ± 10.57 | 53.00 ± 9.78 | 0.299 |
| bili | 3.27 ± 4.54 | 3.08 ± 4.02 | -0.046 |
| albumin | 3.52 ± 0.42 | 3.43 ± 0.43 | -0.222 |
| 事件發生率 | 40.0% | 34.3% | — |
被排除的人平均老了 3.0 歲 （SMD 0.299，超過 0.1 的平衡門檻不少）， 白蛋白也較低（SMD -0.222），而且觀察到的事件率不同。
多重插補在做什麼
多重插補（multiple imputation, MI）的想法可以用三句話講完：
- 用其他變項預測缺失值，並且加上隨機雜訊——雜訊是關鍵，它代表「我們不知道真值」
- 重複 M 次，得到 M 份完整的資料集，各自跑一次分析
- 用 Rubin’s rules 把 M 個結果合併成一個估計與一個標準誤
第 1 步的雜訊常被誤解。如果只填入迴歸預測的條件平均數（單一插補）， 等於宣稱你確定那個值是多少——標準誤會被低估，信賴區間會太窄。
第 3 步的算術是：
- 合併估計 = M 個估計的平均
- 總變異 = 平均組內變異 ＋ (1 + 1/M) × 組間變異
組內變異是「如果資料完整，這個估計有多不確定」；組間變異是「因為要插補，多出來多少不確定」。
library(survival)
data(pbc, package = "survival")
d <- pbc
d$log_bili <- log(d$bili)
d$log_copper <- log(d$copper) # 這一欄有 108 個 NA，本頁處理的就是它們
# 一次插補：迴歸預測 + 隨機殘差
impute_once <- function(data, target, predictors) {
obs <- data[!is.na(data[[target]]), ]
fit <- lm(reformulate(predictors, target), data = obs)
need <- is.na(data[[target]])
# 加雜訊才是插補；只填條件平均數會低估標準誤
data[[target]][need] <- predict(fit, newdata = data[need, ]) +
rnorm(sum(need), 0, sigma(fit))
data
}
M <- 20
est <- se <- numeric(M)
for (m in seq_len(M)) {
dm <- impute_once(d, "log_copper", c("age", "log_bili", "albumin"))
fm <- coxph(Surv(time, status == 2) ~ age + log(bili) + albumin +
log(exp(dm$log_copper)) + factor(stage), data = dm)
est[m] <- summary(fm)$coefficients["log(bili)", "coef"]
se[m] <- summary(fm)$coefficients["log(bili)", "se(coef)"]
}
# Rubin's rules
q_bar <- mean(est) # 合併估計
u_bar <- mean(se^2) # 平均組內變異
b <- var(est) # 組間變異
total <- u_bar + (1 + 1/M) * b # 總變異
fmi <- ((1 + 1/M) * b) / total # 缺失資訊比例
# 實務上一行就好：
# library(mice); imp <- mice(d, m = 20); pool(with(imp, coxph(...)))驗證環境：R 4.6.0 + survival 3.8.6。這裡刻意用 base R 手寫，讓你看到 MI 的每一步；實務上用 mice 套件。
import numpy as np
import statsmodels.api as sm
from statsmodels.imputation.mice import MICEData
pbc = sm.datasets.get_rdataset("pbc", "survival").data
predictors = ["age", "log_bili", "albumin", "log_copper"]
df = pbc.assign(log_bili=lambda x: np.log(x["bili"]),
log_copper=lambda x: np.log(x["copper"]),
status=lambda x: (x["status"] == 2).astype(int))[
predictors + ["time", "status"]]
# MICEData 會做鏈式方程插補
imp = MICEData(df)
est, se = [], []
for _ in range(20):
imp.update_all()
fitted = sm.PHReg(imp.data["time"], imp.data[predictors],
status=imp.data["status"]).fit()
# PHReg 的 params 是 numpy 陣列不是 Series，所以位置要自己查出來——
# 寫死 params[0] 會跟著 predictors 的順序跑，給你的是 age 不是 log(bili)。
# R 端的模型另外校正了 factor(stage)，兩邊接近但不會逐位相同。
k = predictors.index("log_bili")
est.append(fitted.params[k])
se.append(fitted.bse[k])
q_bar = np.mean(est)
u_bar = np.mean(np.square(se))
b = np.var(est, ddof=1)
total = u_bar + (1 + 1/len(est)) * b
fmi = ((1 + 1/len(est)) * b) / totalPython 端用 statsmodels 的 MICEData，或 sklearn 的 IterativeImputer（後者預設只做單一插補，要自己重複並合併）。
結果：兩種做法的比較
| 做法 | 分析人數 | log(bilirubin) 的 HR | 95% CI |
|---|---|---|---|
| 完整個案分析 | 310 | 2.357 | 1.897–2.929 |
| 多重插補（M = 20） | 418 | 2.289 | 1.893–2.767 |
點估計很接近，但信賴區間變窄了——因為插補讓那 108 個人的其他資訊 （年齡、膽紅素、白蛋白都有值）回到了分析裡。
缺失資訊比例（FMI）是 0.052，意思是：這個估計的不確定性裡， 只有 5.2% 來自缺失，其餘來自樣本本身。FMI 高的時候， 插補模型設定得好不好就會嚴重影響結論；FMI 低的時候比較穩健。
還有第三個模型：乾脆不放 copper。 R 腳本其實也跑了它——同一條 Cox，但把 copper 整個拿掉： 412 人、HR 2.459 （2.083–2.903）， 三者之中信賴區間最窄。
但它不能跟上面兩列並排比較。 上面兩列問的是同一個問題（「在 age、bilirubin、albumin、copper、stage 都校正的情況下，bilirubin 的 HR 是多少」），只是處理缺失的方式不同； 不放 copper 的模型換掉了校正變項的集合，因此估計的是另一個量。 它的信賴區間比較窄，不是因為它比較有效率，而是因為它用接受 copper 造成的殘餘干擾換來的—— 把一個共變項刪掉當然可以讓 n 從 310 回到 412， 但你同時也放棄了對它的校正。
這一列在這裡的角色是敏感度分析：如果連換掉校正集合都不會讓結論翻轉， 那結論對「copper 到底要怎麼處理」就相對穩健。它不是第三種缺失值處理方法。
常見誤用
| 誤用 | 為什麼錯 |
|---|---|
| 直接跑模型，讓套件自己丟掉有缺失的列 | 那是完整個案分析，是有假設的選擇，而且不會警告你 |
| 不報告排除了幾人、那些人長什麼樣 | 那張比較表是判斷偏誤方向的唯一依據 |
| 用平均數 / 中位數填補 | 單一插補低估變異，且扭曲變項間的相關結構 |
| 只做一次插補就分析 | 沒有組間變異，標準誤被低估 |
| 插補模型不放結果變項 | 會系統性稀釋暴露與結果的關聯 |
| 對 MNAR 用標準 MI 就當處理完了 | MI 假設 MAR；MNAR 需要敏感度分析（如 delta 調整） |
| 用檢定判斷「是不是 MCAR」就決定做法 | MAR 與 MNAR 在資料上無法分辨，該靠對資料產生過程的理解 |
| 把插補後的資料當成真實觀測值再做其他分析 | 插補值帶著不確定性，必須經過 Rubin’s rules 才有正確的標準誤 |
重跑本頁的所有數字
/opt/homebrew/bin/Rscript figures/scripts/B8-01-missing-data.R讀讀看這張圖
答案取自產生本頁圖表的同一份統計輸出，不是另外打上去的。
這份資料集有 418 列，而 Cox 模型印出來的分析人數是 310。中間那個差額該怎麼讀？
看答案與解析
正確答案: 那 108 是被套件安靜丟掉的——完整個案分析是一個有假設的選擇，不是預設
模型印出來的是留下來的那些列。差額 108 來自 coxph() 遇到任一共變項有缺失就整列跳過，而它不警告、不報錯。106 是拒絕隨機分配、進入登錄組的人數，兩個數字很接近但不是同一群——被丟掉的是「這個模型用到的五個變項至少缺一個」的那些列，只是與登錄組大幅重疊。124 是完整個案這一群裡觀察到的事件數；Cox 模型不會因為沒發生事件就排除任何人，被設限的人照樣貢獻風險集。要判斷丟掉這些列要不要緊，看的是他們與留下來的像不像，不是丟了幾列。
被排除的那 108 列與留下來的 310 列有一張比較表，三個變項各有一個標準化平均差。哪一個說法對？
看答案與解析
正確答案: 年齡的標準化平均差是 0.299，遠超過常用的平衡門檻，被排除的是比較老的一群
年齡那一列的 0.299 已經遠離平衡，被排除的人平均老了三歲左右。膽紅素 -0.046 確實接近零，但一張表要整列一起讀——年齡不平衡、白蛋白不平衡，而且兩群觀察到的事件率不同，用其中一列去替整張表背書會漏掉後面那兩件事。白蛋白 -0.222 那一格則是方向被說反了：被排除的人白蛋白比較低，不是比較高，而符號的意義取決於哪一群當減數，這正是這類表最容易讀錯的地方。完整個案分析不只是損失樣本，它讓世代觀察到的事件率跟著位移；這張表幾乎沒有論文會報，但它是判斷「掉了這些要不要緊」的唯一依據。
多重插補用 Rubin's rules 把不確定性拆成兩塊，這個分析的總變異是 0.0094。哪一個說法對？
看答案與解析
正確答案: 組內變異是 0.008878，它衡量的是資料完整時這個估計本來就有的不確定性，而插補多加的那一塊相對很小
總變異等於平均組內變異，加上一個略大於 1 的因子乘上組間變異。組內變異 0.008878 是「假設資料完整、這個估計本來就有的不確定性」，組間變異 0.000463 是「因為要插補才多出來的那一塊」，兩者差了一個量級以上，所以總變異 0.0094 幾乎全部來自樣本本身。0.051889 不是變異而是缺失資訊比例，它衡量的是插補貢獻的那一塊佔總變異的比例，把它讀成一個變異值，單位就錯了。提醒一句：這一頁的插補是刻意簡化的寫法，會系統性低估組間變異，所以這個比例只能當機制示範，不能當成這份資料的實證結果。
多重插補之後，log(bilirubin) 那個風險比的信賴區間比完整個案分析更窄。為什麼？
看答案與解析
正確答案: 因為插補把分析拉回全部 418 列，被丟掉的那些人在年齡、膽紅素、白蛋白上的資訊回到了模型裡
重點是被丟掉的那些列並不是「什麼都沒有」：他們的年齡、膽紅素、白蛋白都有值，只是 copper 或 stage 缺了一格。完整個案分析為了一格缺失把整列丟掉，插補把那一格補起來、讓整列回到 418 的分析裡，多出來的資訊就是區間變窄的來源。說「進入模型的還是 310」是單一插補的誤解，插補之後跑模型的分析人數確實變了。412 那一列是另一回事：它把 copper 整個拿掉，換掉的是校正變項的集合，因此估計的是另一個量；它的區間比較窄不是因為比較有效率，而是用接受殘餘干擾換來的，不能跟前兩列並排比較。
完整個案分析與多重插補的點估計很接近。有人因此說「缺失集中在一個與暴露無關的變項上，所以無所謂」。這份資料支持這個解釋嗎？
看答案與解析
正確答案: 不支持。copper 與 bilirubin 的 Spearman 相關是 0.63，兩者並不是各走各的，點估計接近是跑出來的結果、不是可以事先推出來的
0.63 是一個中等偏強的單調相關：copper 缺的時候，它與 bilirubin 共享的那一部分資訊也一起缺了，所以「與暴露無關」在這份資料上不成立。膽紅素那一列的標準化平均差確實接近零，但它說的是另一件事——被排除與留下的人在膽紅素的平均上沒有明顯差別，那既不等於缺失機制與暴露無關，也不保證校正之後的係數不會動。0.26 是缺失比例，而比例兩個方向都證明不了：它不能證明沒有偏誤，也不能像那個選項說的那樣保證一定有偏誤——同樣的比例配上不同的缺失機制，結果可以差很多，而這一頁的兩個估計就是接近的。換一個缺失結構——缺失率更高、或缺的是校正之後仍與結果強烈相關的變項——兩種做法就可以差很多。接近與否是結果，不是可以事先假設的前提。
缺失分布圖上方有一整塊實心紅色，橫跨好幾個檢驗欄位。作者說那不是雜訊。哪一個說法對？
看答案與解析
正確答案: 那是 106 位拒絕隨機分配、進入登錄追蹤的病人，試驗檢驗依 protocol 本來就不會做，所以是結構性缺失
那一塊是 106 乘上好幾個檢驗欄位。這群人拒絕隨機分配、進入登錄追蹤，而那些檢驗是試驗 protocol 的一部分——所以缺失的原因是可以寫下來的，這正是判斷它接近 MAR 的依據，而那個判斷來自我們知道資料是怎麼產生的，不是從資料本身看出來的。312 是進入試驗的人數，圖上他們那一段大致是白的，把紅塊讀成試驗組會讓機制的方向整個反過來。血小板則是那張分組表上的例外：它在登錄組也只缺個位數，屬於一般的零星缺失，撐不起圖上那一大塊。缺失有沒有結構、結構是什麼，決定的是後面該用完整個案、插補、還是敏感度分析。
用到這個方法的章節
延伸觀看
Survival Analysis [Simply Explained]素材來源與授權
本頁為原創內容