R 語言統合分析入門 (An Introduction to Meta-Analysis in R)

用 R 進入統合分析的第一扇門

1.1 本章學習目標

讀完本章後,你應該能夠:

  1. 說明統合分析 (meta-analysis) 在臨床醫學與公共衛生研究中的角色。
  2. 分辨系統性回顧 (systematic review) 與統合分析的差異。
  3. 理解效果量 (effect size)、信賴區間 (confidence interval)、權重 (weight) 與森林圖 (forest plot) 的基本概念。
  4. 使用 R 4.5.2 建立一個小型臨床資料集,計算風險比 (risk ratio),並畫出第一張森林圖。

如果你第一次聽到「統合分析」覺得它像某種統計界的火鍋,那你其實已經抓到一半精神了:我們把不同研究的結果放進同一個鍋裡,但不是隨便攪一攪就上桌。統合分析的重點,是用透明、可重複、具統計意義的方法,把多個研究的證據整合成比較穩定的推論。

1.2 為什麼醫學研究需要統合分析?

臨床研究常常不是一錘定音。今天一篇隨機對照試驗 (randomized controlled trial, RCT) 說治療有效,明天另一篇研究可能說效果普通,第三篇研究甚至因為樣本數小而看不出差異。這不是研究者在跟你玩躲貓貓,而是醫學資料本來就會受到樣本數、族群、追蹤時間、照護場域與測量方式影響。

統合分析 (meta-analysis) 是一種定量整合方法,用來合併多個研究的效果估計。它通常建立在系統性回顧 (systematic review) 之上。系統性回顧負責用事先定義好的方法找文獻、篩選研究、評估偏差風險;統合分析則進一步把可合併的數據放進統計模型中,估計整體效果。

用白話說:

  • 系統性回顧問:「我們找到哪些可信的研究?」
  • 統合分析問:「如果把這些研究的數字放在一起,整體證據怎麼說?」

在實證醫學 (evidence-based medicine) 中,這兩者常常是臨床指引、衛生政策與醫療科技評估的重要基礎。不過,統合分析不是「把越多研究加起來就越準」的魔法。若原始研究品質差、族群差異太大、結局定義不一致,合併結果可能只會得到一鍋看似精緻但其實味道很困惑的湯。

1.3 第一個臨床情境:低劑量 aspirin 預防心血管事件

假設我們關心一個臨床問題:

對於高心血管風險成人,低劑量 aspirin 是否能降低 1 年內主要心血管事件?

這是一個簡化過的教學範例。實務上,我們會先定義 PICO:

  • 族群 (population):高心血管風險成人
  • 介入 (intervention):低劑量 aspirin
  • 對照 (comparison):安慰劑或標準照護
  • 結局 (outcome):1 年內主要心血管事件

本章先假設我們已經完成文獻搜尋與篩選,得到 6 個研究。每個研究都報告 aspirin 組與對照組的事件人數與總人數。

Code
aspirin <- tibble::tibble(
  study = c("Taipei Cardio Trial", "Taichung Vascular Study",
            "Kaohsiung Prevention Trial", "Tainan Primary Care Study",
            "Hualien Community Trial", "Chiayi Heart Study"),
  events_tx = c(18, 24, 31, 15, 21, 17),
  total_tx = c(520, 610, 700, 430, 500, 460),
  events_ctrl = c(32, 41, 45, 25, 34, 28),
  total_ctrl = c(515, 600, 690, 425, 490, 455)
)

kable(
  aspirin,
  col.names = c("研究", "Aspirin 組事件數", "Aspirin 組總人數",
                "對照組事件數", "對照組總人數"),
  digits = 3
)
研究 Aspirin 組事件數 Aspirin 組總人數 對照組事件數 對照組總人數
Taipei Cardio Trial 18 520 32 515
Taichung Vascular Study 24 610 41 600
Kaohsiung Prevention Trial 31 700 45 690
Tainan Primary Care Study 15 430 25 425
Hualien Community Trial 21 500 34 490
Chiayi Heart Study 17 460 28 455

1.4 效果量:把每個研究轉成同一種語言

統合分析的第一步,是為每個研究計算效果量 (effect size)。效果量是研究結果的共同貨幣。對二元結局 (binary outcome),常見效果量包括勝算比 (odds ratio)、風險比 (risk ratio) 與風險差 (risk difference)。

本章使用風險比 (risk ratio, RR):

\[ RR = \frac{\text{介入組事件風險}}{\text{對照組事件風險}} \]

RR < 1,表示介入組事件風險較低;若 RR > 1,表示介入組事件風險較高;若 RR = 1,表示兩組風險相同。因為 RR 的抽樣分布通常在對數尺度較接近常態分布 (normal distribution),統合分析常先處理 log(RR),最後再轉回 RR。

Code
aspirin_effects <- aspirin |>
  mutate(
    rr = (events_tx / total_tx) / (events_ctrl / total_ctrl),
    log_rr = log(rr),
    se_log_rr = sqrt(1 / events_tx - 1 / total_tx +
                       1 / events_ctrl - 1 / total_ctrl),
    ci_low = exp(log_rr - 1.96 * se_log_rr),
    ci_high = exp(log_rr + 1.96 * se_log_rr)
  )

kable(
  aspirin_effects |>
    transmute(
      研究 = study,
      `RR` = rr,
      `95% CI 下限` = ci_low,
      `95% CI 上限` = ci_high,
      `log(RR) 標準誤` = se_log_rr
    ),
  digits = 3
)
研究 RR 95% CI 下限 95% CI 上限 log(RR) 標準誤
Taipei Cardio Trial 0.557 0.317 0.980 0.288
Taichung Vascular Study 0.576 0.352 0.941 0.251
Kaohsiung Prevention Trial 0.679 0.435 1.060 0.227
Tainan Primary Care Study 0.593 0.317 1.109 0.319
Hualien Community Trial 0.605 0.356 1.028 0.270
Chiayi Heart Study 0.601 0.333 1.082 0.300

這裡的標準誤 (standard error, SE) 描述估計值的不確定性。標準誤越小,代表該研究對效果估計的資訊量越大。信賴區間 (confidence interval, CI) 則提供一段合理範圍:若研究可以在相同條件下重複很多次,95% 信賴區間的建構程序約有 95% 機率涵蓋真實效果。

請注意,這句話有點繞,但很重要:95% CI 不是「真實值有 95% 機率在這個已觀察區間裡」。在傳統頻率學派 (frequentist) 架構下,真實值被視為固定,隨機的是資料與區間。統計課常在這裡讓大家眉頭打結,沒關係,眉頭打結也是學習曲線的一部分。

1.5 權重:不是每篇研究都同樣有發言權

在統合分析中,樣本數較大、估計較精確的研究通常應該有較大的權重 (weight)。最常見的做法之一是反變異數法 (inverse-variance method):變異數越小,權重越大。

\[ w_i = \frac{1}{SE_i^2} \]

在固定效應模型 (fixed-effect model) 中,我們假設所有研究都在估計同一個真實效果。因此,合併後的 log(RR) 是各研究 log(RR) 的加權平均:

\[ \hat{\theta}_{FE} = \frac{\sum_i w_i \hat{\theta}_i}{\sum_i w_i} \]

其中 \(\hat{\theta}_i\) 是第 \(i\) 個研究的 log(RR)\(w_i\) 是其權重。

Code
aspirin_effects <- aspirin_effects |>
  mutate(
    weight = 1 / se_log_rr^2,
    weight_percent = 100 * weight / sum(weight)
  )

pooled_log_rr <- with(
  aspirin_effects,
  sum(weight * log_rr) / sum(weight)
)

pooled_se <- sqrt(1 / sum(aspirin_effects$weight))

pooled <- tibble::tibble(
  rr = exp(pooled_log_rr),
  ci_low = exp(pooled_log_rr - 1.96 * pooled_se),
  ci_high = exp(pooled_log_rr + 1.96 * pooled_se)
)

kable(
  aspirin_effects |>
    transmute(
      研究 = study,
      `RR` = rr,
      `95% CI` = sprintf("%.2f to %.2f", ci_low, ci_high),
      `權重 (%)` = weight_percent
    ),
  digits = 2
)
研究 RR 95% CI 權重 (%)
Taipei Cardio Trial 0.56 0.32 to 0.98 14.71
Taichung Vascular Study 0.58 0.35 to 0.94 19.44
Kaohsiung Prevention Trial 0.68 0.44 to 1.06 23.64
Tainan Primary Care Study 0.59 0.32 to 1.11 11.96
Hualien Community Trial 0.61 0.36 to 1.03 16.71
Chiayi Heart Study 0.60 0.33 to 1.08 13.53
Code
kable(
  pooled,
  col.names = c("合併 RR", "95% CI 下限", "95% CI 上限"),
  digits = 3
)
合併 RR 95% CI 下限 95% CI 上限
0.606 0.488 0.753

在這個範例中,合併 RR 小於 1,表示 aspirin 組的主要心血管事件風險較低。這不等於我們已經完成臨床建議,因為實務上還要看出血風險、研究品質、異質性 (heterogeneity)、絕對風險差、病人偏好與成本。但作為第一章,這張成績單已經很有精神了。

1.6 森林圖:統合分析的招牌視覺化

森林圖 (forest plot) 是統合分析最具代表性的圖。每一列通常是一個研究,點代表效果估計,橫線代表信賴區間。垂直線通常是無效果線 (line of no effect),對 RR 而言是 1。最下方常放合併估計。

Code
forest_data <- bind_rows(
  aspirin_effects |>
    transmute(
      study, rr, ci_low, ci_high,
      weight_percent,
      estimate_type = "Study"
    ),
  tibble::tibble(
    study = "Pooled fixed-effect estimate",
    rr = pooled$rr,
    ci_low = pooled$ci_low,
    ci_high = pooled$ci_high,
    weight_percent = NA_real_,
    estimate_type = "Pooled"
  )
) |>
  mutate(study = factor(study, levels = rev(study)))

ggplot(forest_data, aes(x = rr, y = study)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey45") +
  geom_errorbar(
    aes(xmin = ci_low, xmax = ci_high, color = estimate_type),
    orientation = "y",
    width = 0.18,
    linewidth = 0.8
  ) +
  geom_point(
    aes(size = if_else(is.na(weight_percent), 14, weight_percent),
        color = estimate_type),
    shape = 18
  ) +
  scale_x_log10(
    breaks = c(0.4, 0.6, 0.8, 1, 1.25),
    limits = c(0.38, 1.35)
  ) +
  scale_color_manual(values = c("Study" = "#356AA0", "Pooled" = "#B23A48")) +
  scale_size_continuous(range = c(2.5, 6), guide = "none") +
  labs(
    x = "Risk ratio (log scale)",
    y = NULL,
    color = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Figure 1.1: 低劑量 aspirin 對 1 年主要心血管事件的統合分析森林圖。

如果你看到森林圖覺得「研究結果大致都在左邊」,這代表大多數研究傾向於支持 aspirin 降低事件風險。不過,請不要只看合併菱形或合併點估計就下結論。森林圖還要看各研究的方向是否一致、信賴區間是否很寬、是否有研究與其他研究明顯不同。統計圖不是裝飾品,它是讓我們放慢速度的提醒。

1.7 R 工作流程:從資料到可重複報告

本書使用 R 4.5.2。你可以在終端機確認版本:

Code
/usr/bin/R --version

在後續章節,我們會逐步使用 metametafornetmeta 等套件。不過,第一章故意先用基本公式與 ggplot2 手作一次,原因很簡單:如果你知道套件背後在做什麼,將來看到輸出表格就不會像看到神秘菜單。

1.7.1 Chapter 1 到 Chapter 9 需要的 R packages

本書 Chapter 1 到 Chapter 9 的範例會用到下列 R packages。建議在開始閱讀前先安裝,之後每章就能直接執行程式碼。

Package 主要使用章節 用途
dplyr 1-9 資料整理與變數轉換
ggplot2 1-9 圖形與森林圖、漏斗圖、SROC 圖
knitr 1-9 Quarto 表格輸出
meta 2, 3, 5, 8 標準統合分析與常見效果量函數
metafor 2, 3, 4, 5, 6, 7, 9 隨機效應模型、異質性、meta-regression、進階模型
metasens 5 發表偏倚與敏感度分析
mice 6 缺失資料的多重插補
netmeta 8 網絡統合分析 (network meta-analysis)
mada 9 診斷準確度統合分析與 Reitsma model
meta4diag 9 Bayesian 診斷準確度統合分析的進階延伸

若你使用的是全新 R 環境,可以一次安裝:

Code
install.packages(c(
  "dplyr", "ggplot2", "knitr",
  "meta", "metafor", "metasens",
  "mice", "netmeta", "mada", "meta4diag"
))

本書寫作與測試時使用 R 4.5.2。若你的套件版本較舊,建議先更新套件再重跑各章腳本;統計結果通常不會大幅改變,但圖形預設值、警告訊息或表格格式有時會因版本不同而略有差異。

建議每個統合分析專案至少包含以下檔案:

  • 原始資料:保留不可覆寫版本。
  • 分析腳本:例如 scripts/chapter1.R
  • 圖表輸出:例如 figs/chapter1_forest.png
  • 可重複報告:例如本章的 Quarto (.qmd) 檔案。

本章對應的 R 腳本儲存在 scripts/chapter1.R。執行方式如下:

Code
/usr/bin/Rscript scripts/chapter1.R

執行後會在 figs/ 產生森林圖 PNG 檔與效果量表格。這種安排有一個好處:報告可以直接執行,腳本也可以單獨重跑。當審稿人、共同作者,或三個月後的自己問「這張圖怎麼來的?」時,你可以優雅地回答:「都在 repo 裡。」這比「我記得好像是用 Excel 拉一拉」更能保護心臟。

1.8 統合分析的基本流程

雖然本章只示範統計合併,但完整統合分析通常包含以下步驟:

  1. 明確定義研究問題,例如使用 PICO。
  2. 制定研究計畫書 (protocol),包含搜尋策略、納入排除標準與分析計畫。
  3. 系統性搜尋資料庫與灰色文獻 (grey literature)。
  4. 雙人獨立篩選文獻與萃取資料 (data extraction)。
  5. 評估偏差風險 (risk of bias)。
  6. 選擇合適效果量與統計模型。
  7. 評估異質性與可能的偏倚。
  8. 報告結果,並進行敏感度分析 (sensitivity analysis)。

這些步驟看起來很多,但不要被清單嚇到。統合分析真正難的地方,不只是 R code,而是研究問題要問得清楚,資料要整理得誠實,模型要符合臨床脈絡。R 是煎匙,不是整間廚房。

1.9 常見新手誤區

第一,看到 p 值 (p-value) 小於 0.05 就過度興奮。統合分析中,效果大小、信賴區間、研究品質與臨床重要性通常比單一 p 值更重要。

第二,把固定效應模型當成「預設正確」。固定效應模型假設所有研究共享同一真實效果,但臨床研究常有不同族群、劑量、追蹤時間與照護環境。下一章我們會正式介紹隨機效應模型 (random-effects model),讓模型承認研究間可能存在真實差異。

第三,只看合併結果,不看單篇研究。合併估計很方便,但它不是免思考通行證。若某篇大型研究與其他研究方向相反,或某些研究品質明顯較差,解讀就需要更細膩。

1.10 小結

本章建立了統合分析的基本直覺:我們先把每個研究轉換成可比較的效果量,再依照估計精確度給予權重,最後產生合併估計與森林圖。你也已經看過一個完整、可重複的 R 小範例。

下一章會進入標準方法的核心:固定效應與隨機效應統合分析。我們會更仔細地談連續型結局、研究間變異,以及為什麼同一批資料在不同模型下可能得到不同的不確定性估計。

1.11 Glossary

中文 English
統合分析 meta-analysis
網絡統合分析 network meta-analysis
系統性回顧 systematic review
隨機對照試驗 randomized controlled trial, RCT
實證醫學 evidence-based medicine
族群 population
介入 intervention
對照 comparison
結局 outcome
效果量 effect size
二元結局 binary outcome
勝算比 odds ratio
風險比 risk ratio
風險差 risk difference
常態分布 normal distribution
標準誤 standard error
信賴區間 confidence interval
頻率學派 frequentist
權重 weight
反變異數法 inverse-variance method
固定效應模型 fixed-effect model
異質性 heterogeneity
森林圖 forest plot
無效果線 line of no effect
研究計畫書 protocol
灰色文獻 grey literature
資料萃取 data extraction
偏差風險 risk of bias
敏感度分析 sensitivity analysis
p 值 p-value
隨機效應模型 random-effects model