傾向分數配對

這一頁在解決什麼問題

世代研究那一章用的是最直接的作法：把干擾因子一起放進 Cox 模型。這個作法在共變項不多、事件夠多的時候很好用，但它有兩個限制。

第一，它把「誰接受治療」與「誰會發生事件」兩件事混在同一個模型裡估。共變項一多，你其實不知道模型是靠哪些人在支撐那個治療係數——可能兩組在某個共變項上根本沒有重疊，模型仍然會照著線性假設外插出一個數字給你，而且不會警告。

第二，事件數決定你能放幾個變項。Cox 模型那頁提過每個變項至少十個事件的經驗法則。這份資料只有 47 個事件，照這個標準最多放四到五個變項——但臨床上值得校正的變項有十幾個。

傾向分數（propensity score）處理的正是這個困境：把十幾個共變項壓成一個數字，於是「校正」這件事的複雜度不再受事件數限制，只受治療組人數限制。

傾向分數是什麼

定義只有一句：在給定共變項的條件下，一個人接受暴露的機率。

實務上用邏輯迴歸估——結果變項是治療，不是臨床結果：

它之所以有用，是因為一個不太直觀的性質：如果一組共變項足以去除干擾，那麼光是把這一個分數配平，就足以讓那一整組共變項在兩組間平衡。 換句話說，12 維的平衡問題被降成一維。這叫平衡分數（balancing score）性質。

這一頁的例子

沿用世代研究那一章的同一個世代：接受根除性前列腺切除術、術中輸血的病人，暴露是輸血用紅血球的儲存時間（最長那一組 vs 短／中），結果是生化復發。用同一個世代是刻意的——同一個臨床問題，用兩種方法回答。

資料集共 316 位病人，扣掉共變項有缺失的 29 位之後，進入分析的是 287 位（暴露組 99、對照組 188），其中 47 位發生生化復發。傾向分數模型放了 12 個共變項。

手把手配一次

下面這段是用 base R 從頭寫的：邏輯迴歸算分數、貪婪最近鄰配對、SMD 從定義算。實務上這一切 MatchIt 一行就做完了（本頁後面有那一行），但只看過那一行的人，其實不知道 caliper 到底擋掉了誰。

library(medicaldata); library(survival)
data(blood_storage, package = "medicaldata")
set.seed(20260822)

d <- blood_storage
d$treat <- as.integer(d$RBC.Age.Group == 3)      # 1 = 最長儲存
covs <- c("Age", "AA", "FamHx", "PVol", "TVol", "T.Stage", "bGS",
          "PreopPSA", "PreopTherapy", "Units", "sGS", "AnyAdjTherapy")
cc <- d[complete.cases(d[, c("TimeToRecurrence", "Recurrence", "treat", covs)]), ]

# ── 1. 傾向分數：結果變項是「治療」，不是臨床結果 ──────────────────
ps_fit <- glm(reformulate(covs, "treat"), data = cc, family = binomial())
cc$ps  <- fitted(ps_fit)
cc$lps <- qlogis(cc$ps)          # 在 logit 尺度上配，兩端才不會擠成一團

# ── 2. 貪婪最近鄰配對，1:1、不放回、caliper = 0.2 個 logit 標準差 ──
cal  <- 0.2 * sd(cc$lps)
tr   <- sample(which(cc$treat == 1))   # 隨機順序：貪婪法對順序敏感
pool <- which(cc$treat == 0)
pairs <- list()
for (i in tr) {
  if (!length(pool)) break
  dist <- abs(cc$lps[pool] - cc$lps[i])
  m <- which.min(dist)
  if (dist[m] > cal) next              # caliper 擋下來：這個人配不到，就不配
  pairs[[length(pairs) + 1]] <- c(i, pool[m])
  pool <- pool[-m]                     # 不放回：用掉的對照不再參與
}
matched <- cc[unlist(pairs), ]
matched$pair <- rep(seq_along(pairs), each = 2)

# ── 3. SMD 從定義算，分母固定用「配對前」的合併標準差 ──────────────
smd_denom <- sapply(covs, function(v) {
  x <- cc[[v]]; t <- cc$treat == 1
  sqrt((var(x[t]) + var(x[!t])) / 2)
})
smd <- function(dat) sapply(covs, function(v) {
  x <- dat[[v]]; t <- dat$treat == 1
  (mean(x[t]) - mean(x[!t])) / smd_denom[[v]]
})
round(cbind(before = smd(cc), after = smd(matched)), 3)

# ── 4. 效果估計，變異數認配對組 ────────────────────────────────────
coxph(Surv(TimeToRecurrence, Recurrence) ~ treat + cluster(pair), data = matched)

import numpy as np, pandas as pd, statsmodels.api as sm
from lifelines import CoxPHFitter

bs = sm.datasets.get_rdataset("blood_storage", "medicaldata").data
bs["treat"] = (bs["RBC.Age.Group"] == 3).astype(int)
covs = ["Age", "AA", "FamHx", "PVol", "TVol", "T.Stage", "bGS",
        "PreopPSA", "PreopTherapy", "Units", "sGS", "AnyAdjTherapy"]
cc = bs[["TimeToRecurrence", "Recurrence", "treat"] + covs].dropna().copy()

ps = sm.Logit(cc["treat"], sm.add_constant(cc[covs])).fit(disp=0).predict()
cc["lps"] = np.log(ps / (1 - ps))

rng = np.random.default_rng(20260822)
cal = 0.2 * cc["lps"].std()
tr = rng.permutation(cc.index[cc.treat == 1])
pool = list(cc.index[cc.treat == 0])
pairs = []
for i in tr:
    if not pool:
        break
    dist = (cc.loc[pool, "lps"] - cc.loc[i, "lps"]).abs()
    j = dist.idxmin()
    if dist[j] > cal:
        continue
    pairs.append((i, j)); pool.remove(j)

idx = [k for p in pairs for k in p]
matched = cc.loc[idx].copy()
matched["pair"] = np.repeat(np.arange(len(pairs)), 2)
CoxPHFitter().fit(matched[["TimeToRecurrence", "Recurrence", "treat", "pair"]],
                  duration_col="TimeToRecurrence", event_col="Recurrence",
                  cluster_col="pair").print_summary()

配出來的結果：91 對、共 182 人，暴露組有 8 位因為 caliper 內找不到夠近的對照而被排除。

重疊：配對之前要先看的那張圖

配對之前應該先問一個比「配得好不好」更前面的問題：兩組的傾向分數有沒有重疊？ 沒有重疊的區域，無論用什麼統計方法都不可能有可比較的對照——那裡的估計只能靠模型外插。

這份資料的傾向分數落在 0.125 到 0.613 之間；兩組共同支撐區（common support，兩組都有人的區間）是 0.233 到 0.609。落在共同支撐區外面的人，配對會自動把他們丟掉——這是配對相對於迴歸校正的一個實質優點：迴歸不會告訴你它在外插，配對會直接把配不到的人列出來。

平衡怎麼看：SMD 不是 p 值

配對做完之後唯一該檢查的東西是平衡。判準是標準化平均差（SMD），慣例上 |SMD| < 0.1 算平衡。SMD 的定義、為什麼它不隨樣本數膨脹、以及 0.1 這條線的來歷，在Table 1 與標準化平均差有完整說明，這裡不重複。

這裡只強調一件事：不要對配對前後做組間 p 值檢定。 配對之後樣本數變小，p 值本來就會變大——就算完全沒有改善平衡，「配對後 p 值不顯著」這句話一樣成立。它衡量的是樣本數，不是平衡。

配對前有 5 個共變項的 |SMD| 超過 0.1，最大的是 0.156；配對後超過門檻的剩下 0 個，最大值降到 0.067，平均值從 0.067 降到 0.029。

雙重穩健：為什麼「再校正一次」不是多此一舉

上一節說「配對後仍不平衡的重要變項，就在配對之後的模型裡再校正一次」。那句話背後有一個名字， 叫雙重穩健（doubly robust）。

配對或加權需要傾向分數模型是對的；配對後的迴歸需要結果模型是對的。兩個都可能錯。 雙重穩健的性質是：兩個模型只要有一個是對的，估計就仍然一致——你不必賭中哪一個， 你買到的是兩張票裡中一張就好。

這也是為什麼「配對之後就不要再放共變項了，不然是重複校正」這種說法是錯的。 配對之後再校正一次，付出的代價是幾個自由度，換到的是一次補救的機會。實務上這是預設做法， 不是例外處理。

配對的三個選擇

配對不是一個動作，是三個獨立的決定。三種設定在同一份資料上的結果：

三個決定各自在買什麼：

• 1:1 還是 1:k。 對照組人多的時候，1:k 可以增加事件數、提高精確度。但第 k 個最近鄰一定比第一個遠，所以k 越大、平衡越差——上表第三列就是這樣：1:2 讓最大 |SMD| 越過了 0.1。
• caliper（配對半徑）。 慣例是 0.2 個 logit 傾向分數標準差（這份資料是 0.067）。它擋掉的是「找不到夠像的人」的那些暴露組個案。上表第一列丟掉 8 位、換到更好的平衡；第二列不設 caliper，一個人都不丟，但最大 |SMD| 比較大。這是偏誤與樣本數之間直接的交換。
• 有無放回（with / without replacement）。 上面的程式碼是不放回：用過的對照從池子裡拿走。放回的作法允許同一位對照被配給多位暴露組，平衡通常更好，但同一個人被重複計算，變異數估計必須另外處理，而且如果少數幾位對照被反覆使用，實際的有效樣本數會遠小於名目人數。

表上最後兩欄是三種設定各自的 HR 與信賴區間（都用配對組穩健變異數）。它們列在這裡不是要你挑一個，而是因為下一段要求的「全部一起報」，長的就是這個樣子：三個點估計從 0.83 到 1.03，信賴區間彼此高度重疊、全部跨過 1。看起來離 1 最遠的是 1:2, caliper 0.2 SD（0.83，0.40–1.74），而它同時是丟掉最多暴露組（23 位）、平衡最差（最大 |SMD| 0.102、1 個變項越線）的那一個——點估計動得最多的那一個設定，也正是樣本被換掉最多的那一個。

敏感度分析的讀法是看這三個數字散得多開，不是從裡面選一個。在這份資料上它們散得不算開，但三個區間都寬到無法排除中等程度的效應，所以散不散開本身也說明不了什麼。

配對後的樣本已經不是原本的母體

這是最容易被忽略、卻直接影響結論怎麼寫的一點。

1:1 配對（尤其是加了 caliper 的）保留的是每一位配得到對照的暴露組個案，加上跟他們最像的那些對照。所以配對後的樣本代表的是「長得像暴露組的那群人」，估到的量是 ATT（average treatment effect on the treated，對已接受治療者的平均效果），不是全體母體的 ATE。

這三者在效果隨個人特徵變化（效果修飾）時數值會不同，而且回答的是不同的臨床問題。ATT 問的是「這些已經在用藥的人，用對了嗎」；ATE 問的是「如果推廣到所有人會怎樣」。逆機率加權那一頁會示範同一份資料在三種估計量下跑出來的數字。

這一頁的模型丟掉了 8 位暴露組個案，所以嚴格說連 ATT 都不是——它是「配得到對照的那部分暴露組」的 ATT。這句話要寫進 limitation，而且要報出被丟掉的人數。

配對後的變異數

配對之後的兩列不是獨立的——它們是因為傾向分數相近才被放在一起的。理論上這代表變異數估計要考慮配對結構，否則信賴區間會失準。作法是在 Cox 模型裡加上 cluster(pair)，改用穩健（sandwich）變異數。

這張表的四列不是同一個估計量，不要拿 HR 的大小來判定哪種方法比較好。「多變項 Cox 校正（全部人）」是全樣本、給定共變項的 conditional HR；「1:1 配對後」則是在配對樣本、以 ATT（average treatment effect on the treated）為目標族群的 HR。兩者的估計對象與目標族群都不同；而且 Cox 的 HR 具有 non-collapsibility——即使完全沒有未校正的干擾，conditional 與 marginal 的數值也不必相等。數字不一樣本身不是「配對改善了估計」的證據。

一行版的 MatchIt

上面手寫的那三十行，實務上就是這一行：

library(MatchIt)
m <- matchit(treat ~ Age + AA + FamHx + PVol + TVol + T.Stage + bGS +
               PreopPSA + PreopTherapy + Units + sGS + AnyAdjTherapy,
             data = cc, method = "nearest", distance = "glm",
             caliper = 0.2, std.caliper = TRUE, replace = FALSE)
summary(m)                 # 平衡表，含配對前後 SMD
md <- match.data(m)        # 配對後的資料，帶 subclass 與 weights 欄
coxph(Surv(TimeToRecurrence, Recurrence) ~ treat + cluster(subclass), data = md)

MatchIt 4.7.2 配出 92 對、HR 1.04（0.51–2.13），與手寫版的 91 對、HR 1.02 只差在貪婪配對的隨機順序不同。兩邊數字不完全相同是預期內的——貪婪配對對「先處理誰」敏感，這也是為什麼配對腳本一定要 set.seed()。

實務上請用套件：它處理好了 optimal matching、full matching、放回的權重、以及一整套平衡診斷。手寫版的用途是讓你知道套件在做什麼，以及在它給你奇怪結果的時候知道要去看哪裡。

結果怎麼寫

三種模型的點估計都落在 1 附近，信賴區間全部跨過 1。正確的寫法是：在這個世代中，未偵測到紅血球儲存時間與生化復發的關聯。

不可以寫成「儲存時間與復發無關」，也不可以把配對後的 HR 1.02 講成「風險增加」——信賴區間 0.51–2.03 同時容納了實質的保護與實質的危害，這份資料的事件數（配對後 32 例）不足以排除中等程度的效應。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B6-02-psm.R

延伸觀看

▶資料小探 – 傾向分數配對法

▶Propensity scores: Everything you need to know in 5min

▶How Propensity Scores Work | NEJM Evidence

▶Propensity score matching: an introduction

▶Outcome research: Causal inference & Propensity score II 傾向分數

回報這一頁的問題