效果量與變異數

統合分析合併的不是「結果」，是一對數字

一篇統合分析看起來像是把好幾篇研究的結論加起來，但實際被送進模型的東西只有兩個：效果量（effect size / effect measure，每篇研究的效應估計）與它的變異數（variance，這個估計有多不精確）。其餘的一切——樣本數、追蹤長度、事件數——都只透過變異數影響結果。

這件事有兩個直接後果，本頁其餘篇幅都在講它們：

1. 選哪一個效果量是一個真正的決定，不是格式問題。同一批四格表用 RR、OR 或 RD 表達，合併值不同、異質性也不同。
2. 變異數就是權重。研究之間的相對影響力完全由變異數決定，而變異數又由事件數（不是樣本數）主導。

六種常見效果量，性質不一樣

為什麼比值型一定要在 log 尺度上合併

這不是慣例，是抽樣分布的形狀決定的。

比值有一個結構性的不對稱：往下最多只能到 0，往上沒有上限。「風險減半」是 0.50，「風險加倍」是 2.00，兩者在臨床上是對稱的兩件事，但在數線上一個離 1 只有 0.50、另一個離 1 有 1.00。取對數之後兩者對 0 就對稱了。

實際模擬一次會更清楚：每組 100 人、對照組事件率 0.20、治療組 0.10（真實 RR = 0.50），重複 40000 次試驗。

差別不只是好不好看。用同一批模擬資料，各算一個 95% 的 Wald 信賴區間：

• 直接在原尺度上做（RR ± 1.96 × SE）：實際涵蓋真值的比例是 92.9%
• 在 log 尺度上做再取回指數（exp(log RR ± 1.96 × SE)）：95.7%

前者比名目上的 95% 低了約 2.1 個百分點，而且它的下界可以掉到 0 以下——一個「負的相對風險」。所以套件的內部運算全部在 log 尺度上進行，escalc() 給你的 yi 是 log RR 不是 RR，只有最後畫圖與報告時才 exp() 回來。

變異數從哪裡來，又怎麼變成權重

以 log RR 為例，變異數的公式是

其中 $a$、$c$ 是兩組的事件數，$n_1$、$n_2$ 是兩組人數。注意主導項是事件數的倒數，不是樣本數的倒數。 一個收了三萬人但只發生二十件事件的研究，精確度接近一個收了兩百人、也發生二十件事件的研究。統合分析裡「大研究」指的永遠是事件多的研究。

固定效應下的權重就是變異數的倒數，正規化成百分比：

用 metadat::dat.bcg 這 13 篇 BCG 疫苗預防結核的對照試驗實際算一次（其中只有一部分是隨機分配，分配方式在異質性那一頁被當成亞組變項）：

最大的一篇（TPT Madras, 1980，176782 人）光是它一篇就拿走 41.4% 的權重，權重最小的一篇（Comstock & Webster, 1969）只有 0.31%——相差超過一百倍。要注意「權重最小」不等於「人數最少」：後面這一篇收了 4839 人，在這 13 篇裡人數排在中間，比它人數更少的那幾篇權重全都比它高。權重低的原因是它的事件數少，不是它的人數少。這正是森林圖上方塊大小要盯著看的理由，也是下一頁換成隨機效應之後會被大幅改寫的東西。

連續型結果：MD 與 SMD 不能混用

當各研究用同一個工具測同一件事（住院天數、收縮壓、體重），直接合併平均差即可，單位保留，讀者一看就懂。

各研究用不同量表測同一個構念（三種不同的憂鬱量表、四種不同的疼痛評分）時，就得先除掉單位。SMD 的做法是把每篇研究的平均差除以它自己的組內標準差：

用 metadat::dat.normand1999（9 篇中風照護單位的試驗，結果變項是住院天數）同時跑兩種：

兩個都未達統計顯著（信賴區間都跨過 0），但要注意的不是那個，而是 I² 從 99.0% 變成 95.5%。這兩個 I² 不可以直接比較。 I² 是個比例，換效果量的同時也換掉了它的分子與分母：τ² 的單位從「天的平方」變成無單位，各研究的組內變異數也跟著換算方式改變（SMD 的變異數還多含一項估標準差的不確定性）。所以這裡看到的下降不能讀成「研究之間的真實差異變小了」——要回答那個問題，得在同一個尺度上看絕對量（見異質性那一頁）。

零格：一個沒有標準答案的問題

事件很罕見時，某些研究會出現某一格是 0。這時 log RR 或 log OR 是無限大，變異數也是無限大，模型直接卡住。

metadat::dat.nielweise2007 是個乾淨的例子：18 篇比較抗菌塗層中心靜脈導管與一般導管的試驗，結果是導管相關血流感染。其中 6 篇有零格，而這 6 篇裡有 1 篇是兩組都掛零。四種常見處理方式，跑出來是這樣：

這張表有兩件事一起在變，讀的時候要拆開。 前兩列是隨機效應模型（REML），差別只在零格怎麼處理；後兩列的 Mantel-Haenszel 與 Peto 本身只有固定效應版本，所以它們與前兩列之間同時差了「零格處理」與「模型」兩件事。要做「只差零格處理」的對照，看的是前兩列。

同一批資料，合併 OR 從 0.299 到 0.439，差距超過四成。四種都是文獻上看得到的做法，沒有一種是錯的，但論文如果不寫他們用了哪一種，這個數字就無法被複製。

換一個效果量，異質性就換一張臉

最後一個容易被忽略的後果：效果量的選擇會改變研究之間看起來有多不一致。 同一批 BCG 試驗，三種效果量：

理由是這樣：這些試驗的基線結核發生率相差好幾個數量級（都市貧民區與鄉村完全不同）。在比值尺度上，基線率被約分掉了，所以只要疫苗的相對保護力接近，各研究就看起來一致。在差值尺度上，基線率留在效果量裡，基線率高的地方風險差自然大得多，於是研究之間差得很開，而且因為 RD 的變異數很小，那個差異在統計上顯得極其確定。

動手跑一次

library(metafor)
library(metadat)

data(dat.bcg, package = "metadat")

# escalc() 是把原始資料轉成 (yi, vi) 的那一步。
# 回傳的 yi 是 log RR，不是 RR —— 合併全程在 log 尺度上做。
rr <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg,
             data = dat.bcg, slab = paste(author, year, sep = ", "))
head(rr[, c("yi", "vi")])

# 同一批四格表換成 OR 與 RD，看合併值與異質性怎麼變
or <- escalc(measure = "OR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg)
rd <- escalc(measure = "RD", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg)
sapply(list(RR = rr, OR = or, RD = rd), function(d) rma(yi, vi, data = d)$I2)

# 權重就是變異數的倒數（固定效應）
fe <- rma(yi, vi, data = rr, method = "FE")
round(weights(fe), 1)

# ── 連續型結果：MD 與 SMD ──
data(dat.normand1999, package = "metadat")
escalc(measure = "MD",  m1i = m1i, sd1i = sd1i, n1i = n1i,
       m2i = m2i, sd2i = sd2i, n2i = n2i, data = dat.normand1999)
escalc(measure = "SMD", m1i = m1i, sd1i = sd1i, n1i = n1i,
       m2i = m2i, sd2i = sd2i, n2i = n2i, data = dat.normand1999)  # 已含小樣本校正

# ── 零格：四種處理方式 ──
data(dat.nielweise2007, package = "metadat")
d <- dat.nielweise2007

# (1) 預設：只對有零格的研究加 0.5
rma(escalc("OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = d,
           add = 1/2, to = "only0"))
# (2) 不校正 —— 有零格的研究會變成 NA 被丟掉，這正是問題所在
rma(escalc("OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = d, add = 0))
# (3) Mantel-Haenszel：不需要校正
rma.mh(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = d)
# (4) Peto：事件罕見、兩組人數接近時的選擇
rma.peto(ai = ai, n1i = n1i, ci = ci, n2i = n2i, data = d)

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])

# statsmodels 沒有 escalc()，四格表要自己轉成 log RR 與變異數
a, n1, c, n2 = tpos, tpos + tneg, cpos, cpos + cneg      # 各研究的四格
log_rr = np.log((a / n1) / (c / n2))
var    = 1 / a - 1 / n1 + 1 / c - 1 / n2                 # 主導項是事件數

res = combine_effects(log_rr, var, method_re="dl")
print(res.summary_frame())          # 輸出仍是 log 尺度，要自己 exp() 回來
print("pooled RR =", np.exp(res.mean_effect_re))

# 零格要自己處理，否則 log() 會給 -inf
zero = (a == 0) | (c == 0)
a_c, c_c = a.astype(float).copy(), c.astype(float).copy()
a_c[zero] += 0.5; c_c[zero] += 0.5   # 只是最粗糙的一種做法，見上文的三個問題

怎麼讀報表

看到一篇統合分析，這一頁對應 Methods 裡「Statistical analysis」的頭幾句，重點看四件事：

1. 合併的是哪一個效果量，有沒有說為什麼。 只寫「以隨機效應模型合併」而不說合併的是 RR 還是 OR 還是 RD，是實質的漏寫。理由該是機轉或臨床上的，不該是「因為異質性比較低」。
2. 原始資料是怎麼轉成效果量的。 論文只給合併值、不給各研究的 $n$ 與事件數，讀者就無法檢查任何一格。附錄裡該有一張逐篇的四格表或效果量表。
3. 有沒有零格，怎麼處理的。 罕見事件的統合分析而完全沒提零格，通常代表用了套件預設而不自知。
4. OR 有沒有被當成 RR 唸。 摘要寫「風險降低 X%」但表格裡是 OR，是最常見的過度陳述，尤其在事件不罕見的時候。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B7-01-effect-size.R

延伸觀看

▶【Biostatistics Corner】Meta-analysis 基礎篇

▶Understand What a Meta-Analysis is in Less Than 5 Minutes

▶An Introduction to Systematic Review and Meta-analysis

回報這一頁的問題