敏感度與影響分析

兩個不同的問題：敏感度分析與影響分析

審稿人回覆裡出現頻率最高的方法學要求之一是「請補做敏感度分析」。這句話底下其實有兩個不同的問題，而它們需要的圖不一樣：

前者是敏感度分析（sensitivity analysis），它的通論在敏感度分析那一頁；後者是影響分析（influence analysis）。本頁講第二種，因為它是全站唯一會逐篇檢查納入研究的地方；leave-one-out 則同時屬於兩邊——它既是「拿掉一個資料點」，也是「換一個納入條件」。

本頁用的是 metadat::dat.bcg，13 篇 BCG 疫苗預防結核的對照試驗，與森林圖與漏斗圖是同一份資料、同一個合併分析：隨機效應（REML），合併 RR 0.489（95% CI 0.344–0.696，p < 0.001），I² 92.2%，τ² 0.3132，Q = 152.2（df = 12，p < 0.001）。產生本頁數字的腳本用一行斷言把這個合併值釘死，兩頁的 RR 不可能對不起來。

Leave-one-out：一次拿掉一篇

最直觀的一個：把 13 篇裡的每一篇輪流拿掉，各重跑一次隨機效應模型，看合併估計走到哪裡。

三件事值得讀出來：

1. 合併值移動的範圍是 RR 0.4517 到 0.5334，相對於全體的 0.4894，最大的一次相對變動是 9.0%（拿掉 Hart & Sutherland, 1977）。
2. 全部 13 個 leave-one-out 區間都仍然排除 1。 沒有任何單獨一篇研究撐著這個結論。
3. 異質性不會因為刪掉任何一篇而消失：I² 最低是拿掉 TPT Madras, 1980 的 87.0%，最高是拿掉 Aronson, 1948 的 93.2%，全體是 92.2%。移動的幅度遠小於它本身的高度，換句話說異質性不是某一篇造成的——這是異質性那頁的重點在這裡的具體版本。

另外一個常見的技術陷阱：leave1out() 的輸出全部在 log 尺度上，沒有做任何轉換。只把點估計 exp() 回去、忘了區間上下界的話，會得到一張「區間不包含自己的點估計」的表——而那張表看起來不會怪到讓人停下來。

逐研究的影響診斷量

leave-one-out 只看合併值移動多少。迴歸診斷那一套（influence()）把同一件事拆得更細，每一篇研究給一組數字：

離群與影響力是兩件事，這是整組診斷量最要緊的區別。 一篇又小又偏的研究可以是明顯的離群值卻幾乎不影響合併值（權重太小）；一篇巨大而中規中矩的試驗可以完全不離群卻主導整個結果（權重太大）。Cook’s D 之所以是常用的總結量，就是因為它同時吃進這兩者。

標記了零篇，卻有一篇讓合併值明顯移動

metafor 的 influence() 會在輸出最右邊印一欄星號，標出它認為有影響力的研究。它用的規則是四條的聯集，任一條成立就標記：

在這 13 篇上，被標記的研究有 0 篇——一篇都沒有。

而且這不是門檻被調鬆了才沒響：最大的 Cook’s D 只走到門檻的 50.4%，最大的絕對 DFFITS 是 60.1%，最大的 hat 是 44.2%。離最近的那條線都還有一段距離。

Baujat plot：一張圖上的兩個問題

Baujat plot 把上面那組診斷量壓成一張散布圖，兩個軸各代表一個問題：

• 橫軸 = 該研究對異質性的貢獻，也就是標準化殘差的平方 resid²/(τ² + vᵢ)。往右代表它與其他篇不一致。這 13 個橫軸值加起來是 12.1。
• 縱軸 = 拿掉它之後合併值移動的平方，除以 leave-one-out 的變異數。往上代表拿掉它，合併值會動。

四個象限各自的意思：

在這份資料上，落在右上象限的是 6 篇：Ferguson & Simes, 1949、Hart & Sutherland, 1977、Vandiviere et al, 1973、TPT Madras, 1980、Rosenthal et al, 1961、Comstock et al, 1976。左上只有 Frimodt-Moller et al, 1973 一篇，右下只有 Comstock & Webster, 1969 一篇，其餘 5 篇落在左下。

影響最大的那一篇：它不是壞掉，是坐在 moderator 的一端

三個工具指向同一篇：Hart & Sutherland, 1977。Cook’s D 最大（0.229），leave-one-out 移動最大（9.0%），Baujat plot 上它在右上象限的最上方。

先看它是一篇什麼樣的研究——以下每一格都直接來自 dat.bcg，沒有任何額外假設：

這不是一篇有問題的研究。 它是隨機分配的、規模在前段、事件數充足、區間很窄。它之所以影響大，是因為它的效果比合併值強得多（自身 RR 0.237 對上合併的 0.489），而且它有足夠的權重讓這個差別被聽見。

那它為什麼效果特別強？異質性那頁已經給了答案：緯度。BCG 在高緯度地區的效果明顯較好，一般認為與低緯度環境中非結核分枝桿菌的暴露造成部分交叉免疫有關。而 Hart & Sutherland, 1977 位在 52 度——它是全體裡距赤道第 2 遠的一篇，也是緯度四十度以上規模最大的一篇。

把緯度放進模型，這件事會立刻顯現出來：

最後那一列是整頁的重點。在只有截距的模型裡，Hart & Sutherland, 1977 的 studentized residual 是 -1.451——看起來偏離。把緯度放進模型之後，它縮到 -0.664。 換句話說，一旦模型知道這篇試驗是在 52 度做的，它觀察到的 log RR（-1.442）跟模型預測的（-1.262）就沒有多不一樣了。

它從來不是離群值，它只是坐在 moderator 的一端。

發現一篇影響力大之後，該做什麼

一個可以照著走的順序：

1. 先確認不是資料錯誤。 事件數、樣本數、方向有沒有抄反（治療組與對照組互換是最常見的一種），效果量的尺度對不對。這一步是唯一一個「發現問題就直接改」的步驟，而且改的是資料不是納入條件。
2. 看它在哪一個維度上與其他篇不同。 人群、介入強度、對照組拿什麼、追蹤長度、結果如何定義、年代、地點、偏誤風險。這一步靠的是讀原文，不是統計。
3. 把那個維度放進模型。 有足夠篇數就做 meta-regression 或事先指定的亞組分析（每個共變項至少十篇，見異質性那頁）；篇數不夠就寫進討論，說明異質性的可能來源。
4. 只有在有獨立於結果的理由時才排除。 「它偏誤風險高」是理由，「它讓合併值移動」不是。而且排除條件要在計畫書裡事先指定；事後才決定的排除，一定要同時報出納入與排除兩個結果，並說明哪一個是主要分析。
5. 報告時把過程寫出來。 「敏感度分析顯示結果穩健」這一句本身沒有資訊。要寫的是：做了哪些擾動、每一個擾動後的估計是多少、有沒有哪一個讓結論改變。

怎麼讀報表

1. 有沒有做影響分析。 沒有的話，讀者無法知道那個合併值是不是由一兩篇撐起來的。森林圖上的方塊大小是一個粗略的替代品，但它只講權重，不講移動。
2. leave-one-out 有沒有報區間，還是只報點估計。 只給一排點估計，看不出結論會不會翻。
3. 有沒有把「被標記」當成唯一的判準。 一篇都沒被標記卻有一篇讓估計動了將近一成，就是這一頁的資料本身。
4. 有沒有排除研究，以及排除的理由是不是獨立於結果。 事先指定與事後決定，在方法學上是兩件不同的事。
5. 排除之後的結果是主要分析還是敏感度分析。 這兩個位置的意思差很多。
6. 有沒有解釋。 找到影響最大的一篇卻只寫「排除後結論不變」，是把最有資訊量的一段觀察丟掉了。

動手跑一次

library(metafor)
library(metadat)

data(dat.bcg, package = "metadat")
d <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg,
            data = dat.bcg, slab = paste(author, year, sep = ", "))
res <- rma(yi, vi, data = d, method = "REML")

# ── 1. Leave-one-out ───────────────────────────────────────────────────────
# 輸出全部在 log 尺度上：CI 上下界也要 exp()，只 exp() 點估計會得到
# 一張「區間不包含自己的點估計」的表，而那張表看起來不會怪。
l1 <- leave1out(res)
data.frame(omitted = res$slab,
           rr = exp(l1$estimate),
           lb = exp(l1$ci.lb), ub = exp(l1$ci.ub),
           I2 = l1$I2)

# 相對於全體合併值的變動百分比 —— 值得看的就是這一欄
100 * (exp(l1$estimate) - exp(as.numeric(res$beta))) / exp(as.numeric(res$beta))

# ── 2. 逐研究影響診斷量 ────────────────────────────────────────────────────
infl <- influence(res)
infl                 # rstudent, dffits, cook.d, cov.r, tau2.del, QE.del, hat, weight
plot(infl)           # 每個統計量一張小圖，門檻畫成虛線
infl$is.infl         # 星號那一欄的邏輯向量；本例全為 FALSE

# 門檻不是魔法數字，可以自己算：p 是模型參數個數，m 是 Cook's D 的自由度
p <- infl$p; k <- res$k
3 * sqrt(p / (k - p))       # DFFITS 的門檻
qchisq(0.5, df = infl$m)    # Cook's D 的門檻
3 * p / k                   # hat 的門檻

# ── 3. Baujat plot ─────────────────────────────────────────────────────────
# 它一定會畫圖，同時把座標 invisible 回傳。想只要座標就先 pdf(NULL)。
bj <- baujat(res)
head(bj[order(-bj$y), ])    # 縱軸最高的幾篇 = 拿掉會讓合併值動最多的幾篇

# ── 4. 影響最大那一篇，為什麼影響大 ────────────────────────────────────────
top <- which.max(infl$inf$cook.d)
dat.bcg[top, ]              # 它的共變項：緯度、分配方式、四格表

# 把 moderator 放進模型，它的殘差就縮回來了
mr <- rma(yi, vi, mods = ~ ablat, data = d, method = "REML")
infl$inf$rstudent[top]              # 只有截距的模型
influence(mr)$inf$rstudent[top]     # 模型含緯度

# 拿掉它，緯度的斜率與 R^2 各變多少 —— 這才是「刪掉它」的真實代價
rma(yi, vi, mods = ~ ablat, data = d, method = "REML", subset = -top)

import numpy as np
from statsmodels.stats.meta_analysis import combine_effects

# metadat::dat.bcg 的四格表，13 篇。帶在這裡，這一段才自己跑得動。
tpos = np.array([4, 6, 3, 62, 33, 180, 8, 505, 29, 17, 186, 5, 27])
tneg = np.array([119, 300, 228, 13536, 5036, 1361, 2537, 87886, 7470, 1699,
                 50448, 2493, 16886])
cpos = np.array([11, 29, 11, 248, 47, 372, 10, 499, 45, 65, 141, 3, 29])
cneg = np.array([128, 274, 209, 12619, 5761, 1079, 619, 87892, 7232, 1600,
                 27197, 2338, 17825])
log_rr = np.log((tpos / (tpos + tneg)) / (cpos / (cpos + cneg)))
var    = 1 / tpos - 1 / (tpos + tneg) + 1 / cpos - 1 / (cpos + cneg)

# log_rr 與 var 在上面算好了（escalc 的等價步驟）
k = len(log_rr)
full = combine_effects(log_rr, var, method_re="dl")
mu_full = full.mean_effect_re

# leave-one-out 得自己寫迴圈；沒有現成函式
for i in range(k):
    keep = np.arange(k) != i
    r = combine_effects(log_rr[keep], var[keep], method_re="dl")
    pct = 100 * (np.exp(r.mean_effect_re) - np.exp(mu_full)) / np.exp(mu_full)
    print(i, np.exp(r.mean_effect_re), round(pct, 2), r.i2)

# Cook's D、DFBETAS、hat、Baujat 座標在 statsmodels 裡沒有對應物件。
# 硬要做的話得自己實作 metafor 的公式，包括每次重擬合的 tau^2；
# 這一頁的內容實務上請用 R。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B7-07-influence-diagnostics.R

延伸觀看

▶統合分析的異質性（Heterogeneity）

▶統合分析異質性和小樣本研究偏誤

回報這一頁的問題