敏感度分析

「請補做敏感度分析」到底在要求什麼

敏感度分析是把一個本來就可以合理選擇的分析決定換掉，重跑，看結論會不會翻。

它動的是分析決定。這句話裡的每一個字都在排除別的東西，而那三件別的東西經常被混進來：

最後一列特別容易搞混，所以講清楚：敏感度分析不需要多重比較調整。它問的不是「這裡面有沒有一個是真的」， 而是「同一個結論在不同設定下站不站得住」。對它做調整等於把穩定性的問題當成顯著性的問題來回答。

這一頁的資料，以及它不能拿來說的話

本頁用 survival::rotterdam，2,982 筆乳癌登錄資料， 其中接受荷爾蒙治療的有 339 筆。整份資料沒有任何遺漏值 （0 格），所以「換一種遺漏值處理」這個很常見的敏感度分析軸在這裡不存在， 下面四個軸是別的東西。

十六個都站得住的分析

四個軸，每個軸兩個選項，全部組合：

用完全因子而不是隨手挑一串設定，是因為只有因子設計能把散布歸因到軸上。 挑一串設定只會得到一堆點，而「哪一個決定在主導結果」才是要問的問題。

一句話讀完這張表：這份資料對「有沒有關聯」這個問題沒有給出單一答案， 給出答案的是分析設定。

主導這條曲線的，是一個沒有人會寫進方法段的決定

把每個軸的兩個選項各自平均起來相減，得到那個軸把估計值推動了多少（對數尺度）：

主導的是淋巴結數要不要截頂。在上圖的下層可以直接看出來：截頂那一列的深色圓點幾乎全部集中在左半邊。

這件事值得停一下。論文的方法段會寫終點是什麼、校正了哪些共變項——那是讀者看得到的兩個軸， 而它們在這裡分別排第二與第三。排第一的是一個連方法段都不會提的資料整理決定。

主分析落在中段，而那正是重點

事先指定的主分析是四個軸上各自最少加工的那個選項： 無復發存活、核心共變項、淋巴結數原值、全追蹤。

它的風險比是 0.884（95% CI 0.760–1.028，p = 0.111）， 未達統計顯著，在十六個設定裡排第 12 名。同一條曲線上有 7 個設定是顯著的。

主分析既不是最保守的那一個，也不是最顯著的那一個。它的正當性不來自它在曲線上的位置， 來自它是在看到這條曲線之前就決定好的。 如果正當性來自位置，那麼：

• 挑最顯著的那一個，理由可以寫成「這是效果最明確的設定」
• 挑最保守的那一個，理由可以寫成「我們採取了最謹慎的做法」

兩句話都寫得出來，而且兩句都是在看到結果之後寫的。事後挑一個保守的來假裝謹慎， 與事後挑一個顯著的，是同一個動作。

「敏感度分析顯示結果穩健」是一句空話

這一句在論文裡出現的頻率極高，而它本身不帶任何資訊：沒說做了哪些擾動、 沒說每一個擾動後的估計是多少、沒說有沒有哪一個讓結論改變。

報告時該給的是三件東西：

1. 做了哪些擾動，逐一列出，包括那些沒有改變結論的
2. 每一個擾動後的效果估計與區間，不是只給 p 值、不是只給平衡表、更不是只給一句話
3. 有沒有哪一個讓結論改變，如果有，說明是哪一個、以及為什麼主分析仍然是主分析

穩健不等於對

上面那條曲線是「不穩健」的樣子。但反過來的情形更需要小心：十六個設定全部顯著， 也不能證明結論正確。

理由在本頁的資料裡就看得見。適應症干擾存在於全部十六個設定中， 換結果定義、換共變項、換編碼、換時間窗，沒有任何一個動作碰得到它。 一組整齊一致的曲線，可以整組偏在同一個方向。

影響診斷那一頁用統合分析講過同一件事： 一份合併分析可以在納入研究全都有同一個系統性偏誤時，對逐篇排除完美穩健—— 因為那個偏誤不會因為少了一篇而消失。

敏感度分析能排除的，只有它實際擾動過的那些東西。 沒被擾動的假設，它一個字都沒說。

事後挑最顯著的代價，量出來

前面說「事後挑與事先指定是不同的東西」。那個差別有多大，可以直接量。

做法：把暴露欄位隨機打散，真實關聯因此正好是零，然後重跑同一組 16 個設定。 重複 4,000 次（種子 20260826）， 每一次記兩個數字——最小的那個 p 值，以及事先指定那個設定的 p 值。

• 只報最顯著的那一個：17.9% 的模擬會宣稱顯著
• 只報事先指定的那一個：6.5%

差距是 11.4 個百分點。這不需要任何推論就能解讀——暴露是被打散的， 所以上面每一次「顯著」都是假陽性。

同一批設定、同一份資料、同一個真實關聯（零）。唯一改變的是哪一個結果會被寫進論文。 這就是為什麼「事先指定」不是行政手續。

讀論文時看哪裡

1. 敏感度分析在哪一段。 在附錄一張表裡逐設定列出效果估計，跟在討論裡寫一句「結果穩健」， 是兩件完全不同的事。
2. 有沒有效果估計，還是只有 p 值與平衡表。 只有 p 值的話，看不出結論會不會翻。
3. 主分析與敏感度分析的位置有沒有寫死。 找「pre-specified」「according to the protocol」 這類字眼，以及有沒有預先註冊編號。
4. 列出來的擾動是不是只有沒事的那些。 一組全部通過的敏感度分析， 跟一組被挑過的敏感度分析，在紙上長得一模一樣。
5. 沒有被擾動的是什麼。 這通常比列出來的更重要——見上面「穩健不等於對」。

動手跑一次

library(survival)

d0 <- rotterdam
d0$rfs_time  <- pmin(d0$rtime, d0$dtime)
d0$rfs       <- as.integer(d0$recur == 1 | d0$death == 1)
d0$nodes_cap <- pmin(d0$nodes, 10)
d0$lpgr      <- log(d0$pgr + 1)

grid <- expand.grid(
  outcome  = c("rfs", "os"),
  covars   = c("core", "plus"),
  nodesf   = c("raw", "cap"),
  followup = c("full", "trunc5"),
  stringsAsFactors = FALSE
)

fit_one <- function(d, g) {
  if (g$followup == "trunc5") {
    tt <- 5 * 365.25
    d$rfs <- ifelse(d$rfs_time > tt, 0L, d$rfs); d$rfs_time <- pmin(d$rfs_time, tt)
    d$death <- ifelse(d$dtime > tt, 0L, d$death); d$dtime <- pmin(d$dtime, tt)
  }
  nod <- if (g$nodesf == "raw") "nodes" else "nodes_cap"
  cov <- if (g$covars == "core") c("age", "size", "grade", nod)
         else c("age", "size", "grade", nod, "meno", "lpgr")
  lhs <- if (g$outcome == "rfs") "Surv(rfs_time, rfs)" else "Surv(dtime, death)"
  f <- as.formula(paste(lhs, "~ hormon +", paste(cov, collapse = " + ")))
  summary(coxph(f, data = d))$coefficients["hormon", c("coef", "Pr(>|z|)")]
}

res <- t(sapply(seq_len(nrow(grid)), function(i) fit_one(d0, as.list(grid[i, ]))))
grid$hr <- exp(res[, 1]); grid$p <- res[, 2]
print(grid[order(grid$hr), ], row.names = FALSE)

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

url = "https://vincentarelbundock.github.io/Rdatasets/csv/survival/rotterdam.csv"
d0 = pd.read_csv(url)
d0["rfs_time"] = np.minimum(d0["rtime"], d0["dtime"])
d0["rfs"] = ((d0["recur"] == 1) | (d0["death"] == 1)).astype(int)
d0["nodes_cap"] = d0["nodes"].clip(upper=10)
d0["lpgr"] = np.log(d0["pgr"] + 1)
size = pd.get_dummies(d0["size"], prefix="size", drop_first=True).astype(float)
d0 = pd.concat([d0, size], axis=1)
size_cols = list(size.columns)


def fit_one(outcome, covars, nodesf, followup):
    d = d0.copy()
    if followup == "trunc5":
        tt = 5 * 365.25
        d["rfs"] = np.where(d["rfs_time"] > tt, 0, d["rfs"])
        d["rfs_time"] = d["rfs_time"].clip(upper=tt)
        d["death"] = np.where(d["dtime"] > tt, 0, d["death"])
        d["dtime"] = d["dtime"].clip(upper=tt)
    nod = "nodes" if nodesf == "raw" else "nodes_cap"
    cov = ["age", "grade", nod] + size_cols
    if covars == "plus":
        cov += ["meno", "lpgr"]
    t, e = ("rfs_time", "rfs") if outcome == "rfs" else ("dtime", "death")
    cph = CoxPHFitter().fit(d[[t, e, "hormon"] + cov], duration_col=t, event_col=e)
    return float(np.exp(cph.params_["hormon"])), float(cph.summary.loc["hormon", "p"])


rows = []
for g in itertools.product(["rfs", "os"], ["core", "plus"], ["raw", "cap"],
                           ["full", "trunc5"]):
    hr, p = fit_one(*g)
    rows.append((*g, hr, p))

out = pd.DataFrame(rows, columns=["outcome", "covars", "nodes", "followup", "hr", "p"])
print(out.sort_values("hr").to_string(index=False))

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B8-07-sensitivity-analysis.R

回報這一頁的問題