固定效應與隨機效應

兩個模型問的不是同一個問題

最常聽到的說法是「隨機效應比較保守」。這句話碰巧描述了多數情況下的結果，但它把因果搞反了，而且會導出錯誤的操作（「異質性高就換模型讓區間寬一點」）。真正的差別在於兩者對世界的假設不同，因此估計的對象根本不是同一個東西。

固定效應模型（fixed-effect model） 假設所有納入的研究都在估計同一個真值 $\theta$：

研究之間看起來不一樣，純粹是抽樣誤差。它回答的問題是：「這批研究共同指向的那個效應是多少？」

隨機效應模型（random-effects model） 假設每篇研究有各自的真值 $\theta_i$，而這些真值來自一個分布：

它回答的問題是：「這些真實效應的分布，平均落在哪裡？」 多出來的 $\tau^2$ 就是研究間變異數（between-study variance）。

在同一批資料上，兩者差多少

用 metadat::dat.bcg 這 13 篇 BCG 疫苗試驗（與上一頁、B7-04 同一份資料）跑兩次：

兩件事值得停下來看：

1. 信賴區間的寬度變成原本的 4.4 倍（在 log 尺度上量）。這是 τ² = 0.313 被加進每一篇研究的變異數所致。
2. 點估計本身也移動了很多，從 RR 0.650 移到 0.489。這一點常被忽略——很多人以為換模型只會改變區間寬度。它會改變點估計，而且改變的幅度取決於「大研究與小研究的結果一不一致」。

權重被改寫成什麼樣子

固定效應下權重是 $1 / v_i$；隨機效應下是 $1 / (v_i + \tau^2)$。加上一個對所有研究都相同的常數 $\tau^2$，效果是把權重往平均拉平——原本很精確的研究，它的優勢被 τ² 稀釋掉了。

現在上一節那個「點估計為什麼會移動」有答案了：最大的那篇試驗（TPT Madras, 1980，176782 人）本身的 RR 是 1.012，也就是未偵測到保護效果。 固定效應下光是它一篇就佔 41.4% 的權重，把合併值往 1 拉；隨機效應把它降到 10.2%，合併值就退回其餘研究所指的方向。

τ² 有很多種估法，而且它們不會給同一個答案

隨機效應模型多出來的 τ² 是要估的，不是資料裡直接讀得到的。常見的估計式：

在這份資料上（13 篇、異質性很大）差異不算大，因為研究數還算夠。研究數少的時候差異會變得很重要。

DL 在研究數少時錯在哪

用模擬把它看清楚：真實的 τ² 固定在 BCG 資料估出來的值（0.313），各研究的內部變異數從那 13 篇實際觀察到的變異數中重抽，只改變研究篇數 $k$，每種情境跑 1500 次。

平均值會受少數極端估計拉動，故本段以中位數描述單一統合分析較常遇到的典型估計；零估計率則顯示三法都可能把異質性估成 0。

真值是 0.313。在 $k$ = 5 時，DL 的典型估計只有 0.206——低估了約 34%，而 REML 與 PM 各是 0.241 與 0.244。而且 DL 有 6.8% 的機率直接回報 τ² = 0，也就是說出「這幾篇研究之間沒有異質性」——在一個真實 τ² 明明不小的世界裡。

τ² 被低估的直接後果是信賴區間太窄，因為 $\tau^2$ 就加在每一篇研究的變異數上。這就是下一節要修的東西。

Knapp-Hartung 校正修的是什麼

標準的隨機效應信賴區間用的是常態分位數（1.96），並且把估出來的 $\hat\tau^2$ 當成真值在用。當 $k$ 小、$\hat\tau^2$ 本身很不穩的時候，這兩件事都會讓區間太窄。

Knapp-Hartung 調整（Knapp-Hartung adjustment，也叫 Hartung-Knapp-Sidik-Jonkman，HKSJ）做兩件事：改用 $t$ 分布（自由度 $k - 1$）取代常態，並用一個把 $\hat\tau^2$ 的不確定性納入的變異數估計式。

同一份 BCG 資料，各估計式加不加 KH 的差別：

在 13 篇的情況下加寬約一成，看起來不多。研究數更少時它救的東西才明顯——同一組模擬，看信賴區間實際涵蓋真值的比例：

$k$ = 5 時，名目上 95% 的區間實際只涵蓋真值 85.7%；加上 Knapp-Hartung 之後回到 93.4%。換句話說，一篇只納入五篇研究、用 DL 標準區間的統合分析，它的「95% 信賴區間」實際上大約是一個 86% 區間。

那到底該選哪一個

動手跑一次

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 = ", "))

fe <- rma(yi, vi, data = d, method = "FE")     # 固定效應
re <- rma(yi, vi, data = d, method = "REML")   # 隨機效應（現在的預設）
fe; re

# 權重被改寫成什麼樣子 —— 這張表是本頁第三節的來源
# slab 不會變成 d 的欄位（d$slab 是 NULL），它跟著 weights() 的 names 走。
data.frame(study = names(weights(fe)),
           fixed  = round(weights(fe), 1),
           random = round(weights(re), 1))

# 換 tau^2 估計式，看合併值與區間跟著動
for (m in c("DL", "REML", "PM", "SJ", "ML", "HE")) {
  r <- rma(yi, vi, data = d, method = m)
  cat(sprintf("%-5s tau2=%.4f  RR=%.3f (%.3f-%.3f)\n",
              m, r$tau2, exp(r$beta), exp(r$ci.lb), exp(r$ci.ub)))
}

# Knapp-Hartung：改用 t 分布，並納入 tau^2 的不確定性
rma(yi, vi, data = d, method = "REML", test = "knha")

# 用 meta 套件也一樣（它的輸出格式更接近論文表格）
# library(meta)
# metabin(tpos, tpos + tneg, cpos, cpos + cneg, data = dat.bcg,
#         sm = "RR", method.tau = "REML", hakn = TRUE)

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 與變異數
res = combine_effects(log_rr, var, method_re="dl")
print(res.summary_frame())        # 同時給固定效應與隨機效應兩列

# 權重可以自己算 —— 這正是本頁的重點，公式很短
w_fe = 1 / var
w_re = 1 / (var + res.tau2)
print(np.round(100 * w_fe / w_fe.sum(), 1))
print(np.round(100 * w_re / w_re.sum(), 1))

# Knapp-Hartung 要自己實作：t 分布分位數 + 加權殘差平方和
k  = len(log_rr)
mu = np.sum(w_re * log_rr) / np.sum(w_re)
q  = np.sum(w_re * (log_rr - mu) ** 2) / (k - 1)
se = np.sqrt(q / np.sum(w_re))
from scipy.stats import t
lo, hi = mu - t.ppf(0.975, k - 1) * se, mu + t.ppf(0.975, k - 1) * se
print("KH CI:", np.exp(lo), np.exp(hi))

怎麼讀報表

論文裡與這一頁有關的通常是 Methods 的一兩句加上 forest plot 底下的那一行。要找五件事：

1. 用了哪個模型，理由是什麼。 只寫「random-effects model」而沒說為什麼，是常態但不理想；完全沒寫模型，是實質缺漏。
2. τ² 用哪個估計式估的。 沒寫通常是 DL。舊的 review 尤其如此。
3. 有沒有 Knapp-Hartung。 寫法是 “Hartung-Knapp adjustment” 或 “HKSJ”。研究數少於十篇又沒有它，那個信賴區間要自己在心裡放寬。
4. 模型是不是事後才決定的。 出現「因為 I² 高於某個值，故改採隨機效應」這種句子，就是資料驅動的模型選擇。
5. 權重集中在誰身上。 forest plot 上有一個方塊特別大，就把它單獨拿掉再看一次結論——這件事讀者自己做得到，只要論文有附各研究的數據。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B7-02-fixed-random.R

延伸觀看

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

▶Systematic reviews and meta analysis

▶How to do your first meta-analysis from start to finish

回報這一頁的問題