累積統合分析與資訊量

每來一篇就重算一次，就是同一個問題被檢定很多次

系統性回顧不是一次性的產物。題目一旦成立，新的試驗會持續發表，作者（或 Cochrane 的更新流程）就會把新研究加進去、重跑一次合併模型、再看一次信賴區間有沒有排除無效值。累積統合分析（cumulative meta-analysis）只是把這個實務過程整個攤開來畫：把研究依發表年份排序，每加入一篇就重跑一次，於是一篇研究對應一列，每一列都是「到這篇為止，證據長什麼樣子」。

這張圖同時說了兩件相反的話，而它們都要講：

• 往回看，它讓人看見合併估計其實在很早的某一步就不再改變了，後面那些試驗只是把信賴區間收窄。這是累積統合分析最有名的用途。
• 往前看，它是一連串在同一個問題上做的檢定。每加入一篇就問一次「現在顯著了嗎」，問了幾十次，而每一次都用同一個未經調整的門檻。這跟多重比較那一頁講的是同一件事——只是那裡的軸是 endpoint 與次族群，這裡的軸是日曆時間。

本頁把這兩件事在同一份資料上量出來。

資料：靜脈 streptokinase 那批試驗

metadat::dat.lau1992 是累積統合分析的正典教材資料：33 篇比較靜脈注射 streptokinase（鏈激酶）與對照的隨機試驗，結果是死亡，橫跨 1959 到 1988 年、共 36,974 位受試者。它就是 Lau 等人 1992 年那篇 NEJM 論文背後的資料集，資料本身已依年份排好序。

治療組累積 1,892 件死亡，對照組 2,375 件；對照組的合併風險是 12.88%。全部合併之後，隨機效應（REML）的合併相對風險是 0.795，95% 信賴區間 0.746 到 0.847，p < 0.001，相當於相對風險降低 20.5%。異質性很低：I² 為 2.0%，τ² 為 0.00079，Q 檢定 38.5（df = 32，p = 0.199）未達統計顯著。

答案早就在那裡了

先看整張累積森林圖。每一列都不是單一試驗的結果，而是到那一篇為止的全部證據合併起來的結果。

「合併值停止改變」這句話必須有判準，否則就只是看圖說故事。腳本把判準連同定義一起寫進 stats 檔的 stability.criterion 欄位，內容是：最早的那一步 j，使得其後每一步都同時滿足 (a) 累積合併相對風險落在最終合併值的正負 10% 以內，而且 (b) 95% 信賴區間排除 1。

兩個條件都是必要的。只看 (a)，很早就會滿足然後又被打破；只看 (b)，那正是這份資料裡會來回閃動的東西。而且條件要求的是「其後每一步都成立」，不是「在第 j 步成立」——否則一個運氣好的步驟會被下一篇試驗推翻，判準卻已經宣告穩定了。

照這個判準，答案是：

也就是說，在這個題目上，全部受試者的 86.0% 是在合併估計已經不再實質改變之後才被隨機分派的。這是累積統合分析被拿來講「證據累積的效率」時最常引用的那一類數字，而 Lau 等人 1992 年那篇論文正是這個論證的原型。

顯著性掉回去過一次

比「合併值收斂」更能說明重複檢定問題的，是顯著性本身來回閃動。把每一步的信賴區間有沒有排除 1 排成一列，這份資料的樣子是：

信賴區間首次排除 1 是在第 14 步（Austrian，1977 年，累積 4,084 人）。它撐了三步。到第 17 步（N Ger Collab，1977 年），累積合併值變成 0.844，95% 信賴區間 0.709 到 1.006——又跨回 1 了。要到第 19 步才回來，而且從那裡開始才沒有再掉出去。整段累積過程一共穿越無效線 3 次。

「第一次顯著是哪一年」取決於你用什麼尺度

上一節說「首次排除 1 是在 1977 年」，這句話只在本頁採用的模型與效果量之下成立。換一個組合，年份就會動。同一份資料、同樣的判準（信賴區間排除無效值），四種常見組合給出的答案是：

同樣是相對風險，固定效應與隨機效應之間差了 6 年。原因不神秘：固定效應在早期給小型試驗較大的權重，也不為研究間變異付出代價，所以信賴區間比較窄，比較早排除無效值；隨機效應在早期估到的研究間變異不小（前二十步的 I² 有好幾步超過四成），區間因此較寬。至於相對風險與勝算比的差別，則來自事件常見時兩者本來就不等價。

Required information size：一條線，四個假設

單一試驗有樣本數估計：先講清楚要偵測多大的效果、容忍多少型一與型二錯誤，才算得出要收多少人。Required information size（RIS，所需資訊量）把同一套算式搬到統合分析上：這個題目全部加起來，需要多少人才算「檢定力足夠」？

算式就是最普通的兩比例雙尾樣本數，再把兩組加總：

其中 $p_1$ 是假設的對照組事件率、$p_2$ 是它在假設的相對風險降低之下會變成的值。這頁用的四個輸入是：

代進去得到 RIS = 5,250（每組 2,625）。這批試驗在第 20 步（ISAM，1986 年）累積到 6,935 人時越過這條線，其後還有 13 篇試驗、30,039 位受試者。

那句「換一個數字就分開」不該只是聲明，所以把整條 grid 列出來。只改相對風險降低這一個輸入，其餘三個不動：

從 10% 改到 30%，RIS 從 22,020 掉到 2,218，差了將近十倍；達標的年份從 1988 提前到 1973。這條 grid 才是誠實的呈現方式，單獨一條線不是。 一張圖上畫一條垂直線，看起來像是資料的性質；畫五條，讀者才看得到那是一個旋鈕。

動手跑一次

library(metafor)
library(metadat)

data(dat.lau1992, package = "metadat")
d <- dat.lau1992          # 已依年份排序；累積統合分析吃的就是這個順序

# add = 1/2, to = "only0" 只對有零格的那一篇加 0.5（Baroffio 1986 治療組零事件）。
# 這是 metafor 的預設，寫出來是為了讓這個選擇看得見。
e <- escalc(measure = "RR", ai = ai, n1i = n1i, ci = ci, n2i = n2i,
            data = d, slab = paste(trial, year), add = 1/2, to = "only0")
res <- rma(yi, vi, data = e, method = "REML")

cu <- cumul(res)          # 一步一列：每一列都是「到這篇為止」的合併值
forest(cu, atransf = exp, refline = 0,
       xlab = "Cumulative pooled risk ratio (death)")

cud <- as.data.frame(cu)
cud$rr    <- exp(cud$estimate)
cud$ci_lb <- exp(cud$ci.lb)
cud$ci_ub <- exp(cud$ci.ub)
cud$sig   <- cud$ci_ub < 1        # 信賴區間排除 1（保護的那一側）

# 穿越無效線幾次：sig 這個布林序列變號幾次
sum(diff(as.integer(cud$sig)) != 0)

# 合併值哪一步之後不再改變：兩個條件都要，而且要對其後每一步都成立
tol <- 0.10
ok  <- abs(cud$rr / cud$rr[nrow(cud)] - 1) <= tol & cud$sig
min(which(rev(cumprod(rev(as.integer(ok)))) == 1L))

# 換一個尺度或模型，「第一次顯著」的年份就會動
cumul(rma(yi, vi, data = e, method = "FE"))            # 固定效應、相對風險
cumul(rma(measure = "PETO", ai = ai, n1i = n1i, ci = ci, n2i = n2i,
          data = d, method = "FE"))                    # Lau 原文的 Peto 勝算比

# required information size：兩比例、雙尾、兩組人數相加
ris_total <- function(pc, rrr, alpha, power) {
  p1 <- pc; p2 <- pc * (1 - rrr); pbar <- (p1 + p2) / 2
  n_arm <- (qnorm(1 - alpha / 2) * sqrt(2 * pbar * (1 - pbar)) +
            qnorm(power) * sqrt(p1 * (1 - p1) + p2 * (1 - p2)))^2 / (p1 - p2)^2
  2 * ceiling(n_arm)
}
# 不要只算一個值：把假設的降幅掃過一整排，才看得出它有多敏感
sapply(c(0.10, 0.15, 0.20, 0.25, 0.30),
       function(r) ris_total(0.12, r, 0.05, 0.80))

# 異質性調整版本（本頁報了但沒畫，因為 I^2 太低，兩條線會重疊）
ceiling(ris_total(0.12, 0.20, 0.05, 0.80) / (1 - res$I2 / 100))

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

# metadat::dat.lau1992，33 篇，已依年份排序 —— 累積統合分析吃的就是這個順序。
# 只對含零格的那一篇加 0.5，與 R 端 add = 1/2, to = "only0" 相同。
ai  = np.array([1, 4, 20, 69, 22, 19, 26, 13, 7, 6, 11, 4, 38, 37, 25, 1, 63,
                5, 18, 54, 628, 1, 0, 1, 1, 3, 3, 2, 4, 1, 12, 791, 2])
n1i = np.array([12, 21, 83, 373, 219, 164, 264, 102, 53, 55, 49, 14, 302, 352,
                123, 13, 249, 32, 156, 859, 5860, 28, 29, 19, 21, 49, 35, 107,
                52, 13, 191, 8592, 41])
ci  = np.array([4, 7, 15, 94, 17, 18, 32, 29, 3, 6, 9, 1, 40, 65, 31, 3, 51, 5,
                30, 63, 758, 2, 6, 3, 1, 6, 4, 12, 7, 2, 17, 1029, 5])
n2i = np.array([11, 21, 84, 357, 207, 157, 253, 104, 54, 53, 42, 9, 293, 376,
                107, 11, 234, 26, 159, 882, 5852, 24, 30, 19, 23, 49, 29, 112,
                55, 12, 177, 8595, 25])
a, b = ai.astype(float), (n1i - ai).astype(float)
c, d = ci.astype(float), (n2i - ci).astype(float)
zero = (a == 0) | (b == 0) | (c == 0) | (d == 0)
a[zero] += 0.5; b[zero] += 0.5; c[zero] += 0.5; d[zero] += 0.5
log_rr = np.log((a / (a + b)) / (c / (c + d)))
var    = 1 / a - 1 / (a + b) + 1 / c - 1 / (c + d)

# 累積統合分析要自己迴圈：statsmodels 沒有 cumul() 的等價物。
# log_rr 與 var 在上面從四格表算好了（含零格的連續性校正）。
rows = []
for k in range(1, len(log_rr) + 1):
    res = combine_effects(log_rr[:k], var[:k], method_re="dl")
    mu, se = res.mean_effect_re, res.sd_eff_w_re
    rows.append((np.exp(mu), np.exp(mu - 1.96 * se), np.exp(mu + 1.96 * se)))

rr, lb, ub = map(np.array, zip(*rows))
sig = ub < 1
print("穿越無效線次數 =", int(np.sum(np.diff(sig.astype(int)) != 0)))

# required information size：算式簡單，難的是那四個假設不是資料給的
def ris_total(pc, rrr, alpha=0.05, power=0.80):
    p1, p2 = pc, pc * (1 - rrr)
    pbar = (p1 + p2) / 2
    n_arm = (norm.ppf(1 - alpha / 2) * np.sqrt(2 * pbar * (1 - pbar)) +
             norm.ppf(power) * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) ** 2 / (p1 - p2) ** 2
    return 2 * int(np.ceil(n_arm))

print([ris_total(0.12, r) for r in (0.10, 0.15, 0.20, 0.25, 0.30)])

# 注意：statsmodels 的隨機效應只有 DerSimonian-Laird，與上面 R 的 REML 不同，
# 兩邊的數字不會完全一致。想重現本頁的值請用 R。

怎麼讀報表

1. 累積森林圖的每一列是什麼。 它不是單一試驗的結果，而是「到這篇為止」的合併值。把它當成一般森林圖讀，會以為每一篇試驗都很精確。
2. 排序的依據是什麼。 依年份是最常見的，但也有人依樣本數、依品質、依風險偏差排序。依年份以外的排序不叫「證據何時足夠」，那是敏感度分析，要另外講。
3. 有沒有交代「首次顯著」的門檻。 如果只寫「到某年就達到統計顯著」而沒說是哪個模型、哪個效果量、哪個 p 值門檻，這句話幾乎不可驗證。
4. 有沒有 required information size，以及它的四個輸入寫在哪。 沒有寫出對照組事件率與目標效果的 RIS，等於沒有 RIS——那條線的位置幾乎完全由這兩個數字決定。
5. RIS 有沒有做敏感度分析。 一條線是一個假設，一組線才是一個結論。
6. 有沒有畫 TSA 邊界，畫的是哪一族。 O’Brien-Fleming、Pocock、alpha-spending 收緊的速度不同，同一份資料可以給出不同的「穿越了沒有」。
7. 異質性調整過的 RIS 旁邊有沒有 I²。 沒有 I² 就無法判斷那個調整是實質的還是幾乎沒動。
8. 結論的措辭。 「證據已經足夠」與「在這組假設下，累積資訊量已超過 required information size」是兩句不同的話。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B7-06-cumulative-meta.R

回報這一頁的問題