Forest plot 與 funnel plot

Forest plot：每個元素在編碼什麼

森林圖看起來像一堆點和線，但每個視覺元素都在傳遞特定資訊：

方塊大小是最容易被忽略、卻最有資訊量的元素。 它告訴你這個合併結果實際上是被誰主導的。常見的情況是：一篇超大型試驗的方塊大到佔掉大半權重，其餘十幾篇小研究加起來影響有限——這時候「合併了 15 篇研究」的說法就有點誤導，實質上你在看的接近是那一篇大研究的結果。

library(metafor)
library(metadat)

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

# 從四格表算出各研究的 log risk ratio 與變異數
bcg <- escalc(measure = "RR",
              ai = tpos, bi = tneg, ci = cpos, di = cneg,
              data = dat.bcg, slab = paste(author, year, sep = ", "))

# 隨機效應模型（REML）
res <- rma(yi, vi, data = bcg, method = "REML")
res                                   # 合併估計、I²、tau²、Q 檢定

forest(res, atransf = exp,            # 畫在 RR 尺度上
       at = log(c(0.05, 0.25, 1, 4)),
       header = c("Trial", "Risk ratio [95% CI]"))

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])
rr        = (tpos / (tpos + tneg)) / (cpos / (cpos + cneg))
se_log_rr = np.sqrt(1 / tpos - 1 / (tpos + tneg) + 1 / cpos - 1 / (cpos + cneg))

# 以上面算好的 log RR 與其變異數作為輸入
eff = np.log(rr)            # 各研究的 log risk ratio
var = se_log_rr ** 2        # 對應變異數

res = combine_effects(eff, var, method_re="dl")   # DerSimonian-Laird
print(res.summary_frame())
print("I2 =", res.i2)

合併結果是 RR 0.49（95% CI 0.34–0.70），共納入 13 篇研究。圖上每一個方塊、每一條橫線，都是某一篇研究的 log RR 與它的變異數——這兩個數字怎麼算出來、為什麼決定方塊大小的是變異數，在效果量與它的變異數那頁。

固定效應與隨機效應：問的是不同問題

上面那個菱形出自隨機效應模型，而是哪一種模型畫出來的，會改變它的意思：固定效應的菱形估的是一個共同的真實效果，隨機效應的菱形估的是研究各自真實效果所構成的那個分布的平均。所以讀森林圖的第一步，是先弄清楚菱形是哪一種模型給的——順帶注意隨機效應會把權重分配得比較平均，方塊的大小已經不是固定效應下的那組大小了。

一組試驗該用哪個假設、選擇的代價是區間寬多少、用的是哪一個 τ² 估計式，是固定效應與隨機效應模型那頁的主題，本頁不重講。

I²：它衡量的是比例，不是量

森林圖通常會在菱形下方印出 I²（本例 92.2%），旁邊再附上 Q 檢定（本例 Q = 152.2, df = 12, p = < 0.001）。這兩個數字被引用的次數遠多於被讀懂的次數，而且各有一個標準的誤讀，先擋掉才看得懂這張圖。

I² 的定義是研究間變異佔總變異的百分比——它是個比例，所以 I² 高不代表研究之間的效果差異大。要看絕對大小，該看 τ²（研究間變異數，本例 0.313）或預測區間（prediction interval）。而 Q 檢定在研究篇數少時檢定力很低，所以拒絕不了虛無假設，並不等於各研究一致。

這兩個統計量各自真正的意思、I² 那兩個算出來不一樣的估計式、以及 Q 檢定的檢定力到底有多低的模擬，都在異質性：I²、τ² 與預測區間那頁。

而 BCG 這份資料之所以是經典教材，正因為它的異質性有可解釋的來源：試驗地點的緯度。緯度高的地區 BCG 效果明顯較好（一般認為與環境中非結核分枝桿菌的暴露有關）。這正是異質性的正確處理方式——不是想辦法把它壓下去，而是找出它從哪裡來，用 meta-regression 或亞組分析解釋它。

Funnel plot 與不對稱

漏斗圖把每篇研究的效果量（x 軸）對上它的精確度（y 軸，通常是標準誤，且上下顛倒讓精確的研究在上方）。

邏輯是：如果不存在偏誤，精確的大研究會集中在真值附近（漏斗頸），不精確的小研究會對稱地散在兩側（漏斗口），整體呈倒漏斗形。

Egger 檢定把「用眼睛看」形式化：拿每篇研究的效果估計對它自己的標準誤做迴歸，看斜率是不是偏離零。漏斗對稱時斜率應該是零——估計值不隨研究精確與否而變；斜率偏離零，就表示比較小的研究系統性地落在某一側。（在這個參數化下，截距是另一個有用的量：它是外推到「精確度無限大的研究」所得的合併估計，也就是下面提到的 Egger 校正後估計的來源。Egger 原始寫法是把標準常態離差對精確度迴歸，那邊承擔檢定的才是截距——同一個檢定、不同的代數，這也是不同資料來源對「該讀哪個係數」講法不一的原因。）

這份資料的合併結果是 OR 0.46（95% CI 0.31–0.70），Egger 檢定 p < 0.0001。此資料的 I² 為 73.6%，所以研究間變異是解釋此不對稱時必須考慮的可能來源；I² 是相對變異比例，不可單獨用來判定異質性的絕對臨床大小或不對稱的成因。

同一條迴歸的截距——外推到「精確度無限大的研究」——在這份資料上是 OR 1.105，落在無效值的另一側。這是一個外推值，沒有自己的檢定，只用來說明資料在往哪個方向指。

這是統合分析史上最有名的案例之一：一堆小型試驗一致顯示靜脈鎂對心肌梗塞有顯著好處，合併後看起來效果很大；後來的超大型試驗（ISIS-4，收了超過五萬人）卻未偵測到保護效果——在這份資料裡，它的點估計甚至落在 1 的右側。漏斗圖的不對稱在大型試驗發表之前就已經顯示出警訊。

不對稱檢定不只 Egger 一種

Egger 檢定最常被引用，但它既不是唯一的一種，也不是每種結果變項都該用它。實務上會遇到三種檢定加一種敏感度分析，四者針對的是同一件事——報出來的效果會不會跟著研究的大小走——差別在於「大小」用什麼代表，以及怎麼加權。

下表是它們在這份鎂試驗資料（16 篇）上的結果。

Egger 對二元結果的已知偏誤，就是 Peters 存在的理由

回頭看上面那份小研究效應的清單，第四項是「效果量與變異數的數學相依」。對 OR 與 RR 而言，那不是一句抽象的警語：log OR 的點估計與它的標準誤是從同一組四格表的格子算出來的，一篇碰巧觀察到較大效果的研究，同時也會拿到較大的估計標準誤。於是「效果量對標準誤迴歸」在完全沒有任何偏誤的情況下也會有一個非零的斜率，事件率低、兩組人數不平衡、真實效果大的時候尤其明顯——Egger 檢定的第一型錯誤率因此被推高。

Peters 的做法是把預測變數換成一個不含觀察到的效果的量：1 除以總樣本數。權重也跟著換成由事件數與非事件數的邊際總數建出來的有效樣本數，而不是效果量自己的變異數倒數。同一份資料上，Egger 給 p < 0.0001，Peters 給 p = 0.0014——方向一致，但 Peters 的 p 值大了超過一個數量級，差的那一段有一部分就是上面說的人工相關。

三個檢定在同一份資料上不一致

同樣是這 16 篇試驗，Egger 與 Peters 都指出不對稱，秩相關檢定卻沒有（p = 0.4503）。這不是矛盾，而是檢定力的差別：秩相關檢定放棄了線性假設，換來的代價就是要更多篇研究才偵測得到同樣程度的不對稱。秩相關檢定沒有偵測到不對稱，不能拿來反駁另外兩個檢定——它只是沒有偵測到。

把同樣四種做法搬到本頁前半段那份 BCG 資料上，結果全部翻面：

BCG 這一側，三個檢定都未偵測到漏斗圖不對稱。這句話只能寫到這裡：它不等於「沒有發表偏誤」，也不等於「漏斗是對稱的」。這份資料只有 13 篇，剛好貼著 Cochrane 那條不做不對稱檢定的門檻，而在這個篇數下三個檢定的檢定力本來就低——低到「沒偵測到」幾乎是預設結果，不是資訊。順帶一提，BCG 這份資料的 I² 是 92.2%，異質性極高卻測不出不對稱：異質性與漏斗不對稱是兩件不同的事，一個高不蘊涵另一個高。

Contour-enhanced funnel：缺口落在哪一區

一般的漏斗圖只告訴你哪裡空了，不告訴你空掉的那一塊如果真的有研究，它們會不會顯著。而這一格資訊正是判斷「這個缺口像不像發表偏誤」的關鍵。contour-enhanced funnel 就是把顯著性區域直接塗在底圖上。

等高線以無效值為中心是這張圖的全部意義所在。若把它們畫在合併估計周圍（很多套件的預設），塗出來的形狀長得幾乎一樣，但回答的是完全不同的問題——那時它顯示的是「各研究離合併值多遠」，而不是「各研究自己顯不顯著」。

讀法只有一條，但它把五種機制分成了兩堆：

• 缺口落在不顯著區 —— 與發表偏誤一致。沒被發表的，正是那些做不出顯著結果的小研究。
• 缺口落在顯著區 —— 發表偏誤解釋不通，因為沒有人會把顯著的結果壓在抽屜裡。這時要往異質性、小研究的方法品質、效果量與變異數的數學相依，或單純的偶然去找。

這份資料的分布是：

左側 14 篇裡有 6 篇達到名目顯著；右側只有 2 篇，而且沒有一篇顯著。缺的是「小型、方向偏向對照組、而且不顯著」的試驗——落在白色那一塊。這正是發表偏誤會留下的形狀，也是為什麼這份資料在 ISIS-4 出來之前就該讓人不安。

Trim-and-fill：補上去的研究是演算法生成的

trim-and-fill 的演算法分三步：先修剪掉造成不對稱的那幾篇極端研究，用剩下的重估合併值；再以這個新的合併值為對稱軸，把修剪掉的研究鏡射到另一側補回去；最後用補完的完整資料重算。

在這份資料上，演算法補了 7 篇，全部在右側，研究數從 16 變成 23。合併 OR 從 0.465 移到 0.676，95% CI 0.443–1.031——調整後的區間跨過了無效值。原本看起來很篤定的保護效果，在這個補法之下不再達到統計顯著。

這句話要小心讀。它不是「靜脈鎂沒有效」的證明，理由有三個：

1. 補上去的那 7 個點不是被找回來的研究。 它們沒有病人、沒有事件數、沒有作者，只是演算法照對稱性算出來的座標。真正缺的那些研究（如果存在）是什麼樣子，沒有人知道。
2. 它假設不對稱來自缺漏。 若不對稱其實來自異質性——這份資料的 I² 是 73.6%——那麼補進去的研究是在為一個不存在的缺漏做補償，把估計推向一個沒有依據的方向。
3. 它對估計 τ² 的方法很敏感。 換一個隨機效應估計式，補幾篇、補多遠都會變。

所以 trim-and-fill 的輸出該被讀成一句條件句：如果不對稱真的來自缺漏，合併估計最多會移動這麼多。 那是敏感度分析的一個界，不是校正後的答案。

同樣的做法搬到 BCG 資料上只補了 1 篇，RR 從 0.489 移到 0.518（95% CI 0.365–0.736），區間仍未包含無效值。當漏斗本來就沒有明顯的不對稱，這個演算法能做的事情也就很有限——這是它正常的行為，不是它「確認了沒有偏誤」。

動手跑一次

library(metafor)
library(metadat)

data(dat.egger2001, package = "metadat")
mag <- escalc(measure = "OR", ai = ai, n1i = n1i, ci = ci, n2i = n2i,
              data = dat.egger2001, slab = study)
res <- rma(yi, vi, data = mag, method = "REML")

# 1. Egger：效果量對標準誤。承擔檢定的是斜率；
#    截距是外推到無限精確度的「limit estimate」
regtest(res, model = "lm")

# 2. Begg-Mazumdar：標準化效果量與變異數的 Kendall 秩相關
ranktest(res)

# 3. Peters：預測變數換成 1/總樣本數，權重換成事件與非事件的邊際總數。
#    metafor 的 regtest(predictor = "ninv") 只換了預測變數，權重仍是 1/vi，
#    那不是 Peters，所以這裡照原文的迴歸式自己寫
a  <- dat.egger2001$ai
b  <- dat.egger2001$n1i - a
cc <- dat.egger2001$ci
d  <- dat.egger2001$n2i - cc
ni <- a + b + cc + d
wi <- 1 / (1 / (a + cc) + 1 / (b + d))
summary(lm(mag$yi ~ I(1 / ni), weights = wi))   # 讀第二列：斜率

# 4. Trim-and-fill：不是檢定，是敏感度分析
tf <- trimfill(res)
tf                                    # k0 = 補了幾篇、side = 補在哪一側
predict(tf, transf = exp)             # 補完之後的合併估計與區間

# contour-enhanced funnel。refline = 0 是關鍵：等高線要以虛無值為中心，
# 畫在合併估計周圍的話，同樣的圖回答的是另一個問題
funnel(res, atransf = exp, refline = 0,
       level = c(90, 95, 99),
       shade = c("white", "gray88", "gray72"),
       legend = TRUE)

funnel(tf, atransf = exp)             # 補上去的研究會畫成不同的符號

import numpy as np
import statsmodels.api as sm
from scipy.stats import kendalltau

# metadat::dat.egger2001 的鎂試驗，16 篇。連續性校正**只加在含零格的那一篇**，
# 與 R 端 escalc 的預設 add = 1/2, to = "only0" 相同——每一格都加會讓
# 16 篇裡有 14 篇的 yi 與 R 對不上。
ai  = np.array([1, 9, 2, 1, 10, 1, 1, 0, 6, 1, 2, 5, 4, 90, 4, 2216])
n1i = np.array([40, 135, 200, 48, 150, 59, 25, 22, 76, 27, 89, 23, 130,
                1159, 107, 29011])
ci  = np.array([2, 23, 7, 1, 8, 9, 3, 1, 11, 7, 12, 13, 8, 118, 17, 2103])
n2i = np.array([36, 135, 200, 46, 148, 56, 23, 21, 75, 27, 80, 33, 122,
                1157, 108, 29039])
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
yi = np.log((a * d) / (b * c))
vi = 1 / a + 1 / b + 1 / c + 1 / d

# yi = 各研究的 log OR，vi = 其變異數，都在上面算好了
sei = np.sqrt(vi)

# 1. Egger：加權最小平方，權重 1/vi，檢定斜率
egger = sm.WLS(yi, sm.add_constant(sei), weights=1 / vi).fit()
print(egger.tvalues[1], egger.pvalues[1])

# 2. Begg-Mazumdar：Kendall tau
#    metafor 會先把效果量標準化再算，所以數值與 ranktest() 會略有出入
print(kendalltau(yi, vi))

# 3. Peters：預測變數 1/n，權重用事件與非事件的邊際總數
ni = a + b + c + d
wi = 1 / (1 / (a + c) + 1 / (b + d))
peters = sm.WLS(yi, sm.add_constant(1 / ni), weights=wi).fit()
print(peters.tvalues[1], peters.pvalues[1])

# 4. trim-and-fill 與 contour-enhanced funnel：
#    Python 沒有維護中的實作，這兩步回 R 跑

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B7-04-forest-funnel.R

延伸觀看

▶How to Interpret a Forest Plot

▶醫學統計 EP18 統合分析：加權整合多項研究、看懂森林圖

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

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

回報這一頁的問題