Bland-Altman 分析與一致性界限

這一頁要回答的問題

論文的 Results 常常這樣寫：「新方法與參考方法呈高度相關，顯示兩者具有良好的一致性。」前半句是統計，後半句是願望，而且兩句之間沒有任何推論關係。這一頁用的資料會讓那個相關係數看起來完全站得住腳，然後把同一批人的差值攤開來給你看。

一致性（agreement）問的是一件很具體的事：如果我把這個病人的測量方式從 A 換成 B，我拿到的數字會差多少，而那個差值大到會讓我改變處置嗎。 相關係數回答不了它，因為相關係數量的是「兩個方法把病人排出來的順序有多接近」，不是「同一個病人身上兩個讀數差幾個單位」。

讀完這一頁之後，看到一篇方法比較的論文，你應該問得出四個問題：

1. 偏差有多大、往哪個方向——新方法平均而言比參考方法高還是低，高多少。
2. 一致性界限有多寬——而且要換算成臨床單位：這個寬度等於幾個「臨床上有意義的差」。
3. 界限本身有多不確定——界限是估計出來的，它有自己的信賴區間，而重複測量處理錯了，垮的通常是這一項。
4. 差值會不會隨測量值大小而變——如果會，一組固定的界限就不能套用在整個量程上。

資料：Bland 與 Altman 自己那份血壓資料

本頁用 MethComp::sbp：765 筆讀數，來自 85 位受試者，每位受試者由 3 個方法各測 3 次，單位是 mmHg。這是 Bland 與 Altman 用來示範這套方法的原始教學資料，三個方法的意義如下：

這一頁比的是人工觀察者 J 與半自動機器 S，差值一律定義成 S − J：正值代表機器讀得比人高。可以配成 255 個單次配對，或 85 個「每人三次平均」的配對。

誘導：先看散點圖，再看同一批人的差值

左邊那張圖給出的相關係數是 r = 0.817（95% CI 0.732 到 0.878，p < 0.0001），決定係數 r² = 0.668。用全部 255 個單次配對算則是 0.796。

這個數字沒有任何一篇臨床論文會拒絕稱它為「高度相關」——它的 p 值小到不成問題，區間下界還在 0.73 以上。而它對「這台機器能不能拿來取代人工測量」這個問題，貢獻是零。

這張圖為什麼橫軸放兩者的平均

右邊那張圖的橫軸是兩個方法的平均，不是參考方法。這不是美觀問題：差值與參考方法之間本來就會有人為的負相關（參考方法那一次剛好量偏高時，差值就會被拉低），把差值對參考方法作圖，會憑空生出一條向下的斜率，讓你以為有比例性偏誤。兩者的平均是目前對「真值」最好的近似，它把這個人為關聯攤平掉。

差值圖上的三條線是這樣算的：中間那條是差值的平均，也就是偏差 +15.62 mmHg（95% CI +11.54 到 +19.70）；上下兩條是偏差加上與減去 1.96 倍的差值標準差 20.95 mmHg，得到 -25.4 與 +56.7。跨距 82.1 mmHg。

意思是：這台機器換上去之後，個別病人的讀數可能比人工測量低 25 mmHg，也可能高 57 mmHg，而這還是「95% 的人落在裡面」的範圍。

界限有多寬才算寬：三條把它翻譯成臨床的路

界限是一個統計量，「夠不夠好」不是統計問題。要把它變成可以判斷的東西，有三條路，這份資料三條都走得通。

第一條，換算成臨床上有意義的差。 收縮壓上，10 mmHg 是一般認為足以影響處置的量級。界限跨距等於 8.2 個這樣的差。實際看單次配對的分布也是同一個結論：62.7% 的配對差超過 10 mmHg，23.5% 超過兩倍，而單一配對最大差到 109 mmHg。

第二條，看它會不會翻轉決策。 把 140 mmHg 當成分界，用 J 判定為高的有 24 人，用 S 判定為高的有 36 人。逐人比對，共有 14 位（16.5%）僅僅因為換一個量測方法就跨過這條線——其中 13 位是機器判高而人工判不高，1 位反向。切點怎麼決定、切點附近的誤分類要怎麼看，見切點的選擇。

第三條，找一個外部的法規錨點。 ISO 81060-2 criterion 1 對非侵入式血壓計的要求是：與參考方法的平均差不超過 5 mmHg，差值標準差不超過 8 mmHg。這份資料的觀測值是平均差 15.62 與標準差 20.95 mmHg，兩項都大幅未達。

對照組：同一批人、同一個場合，兩位人工觀察者

這份資料還有第三個方法：另一位人工觀察者 R。把 J 與 R 拿來配對，樣本、場合、統計方法全部不變，只換掉「另一方是誰」：

兩位人工觀察者之間的絕對一致性 ICC 是 0.997，界限跨距 8.9 mmHg——大約是人與機器那一組的 11%，換句話說機器那一組寬了大約 9 倍。

這個對照排除掉一整類辯解。界限之所以寬，不是因為血壓本來就會跳、不是因為受試者不合作、不是因為 85 個人太少、也不是因為 Bland-Altman 這個方法對誰都嚴苛。同一批人、同一個下午、同一組公式，兩個人拿同一台水銀血壓計互相對照，差值標準差只有 2.26 mmHg。寬的是這台機器。

重複測量：垮掉的是精確度，不是界限

每位受試者被各測了 3 次，所以「有幾個觀測值」這個問題有三種答法，而三種答法給出三組不同的界限：

單位都是 mmHg。C 是本頁採用的做法。三橫列擺在一起，結論與多數人的直覺不一樣：

把 255 個配對全部當成獨立觀測，界限本身幾乎沒有動。 跨距 79.84 對上正確的 82.12，只窄了 2.77%，換算成絕對值是 2.28 mmHg。畫出來的圖，兩條虛線的位置肉眼分不出來。

垮掉的是界限的精確度。 界限的 95% 信賴區間從正確的上下各 6.60 mmHg，縮成上下各 4.30 mmHg，窄了 34.9%。原因很直接：有效樣本數被灌水 3 倍——資料裡只有 85 個獨立的人，卻按 255 個獨立觀測去算標準誤。

真正嚴重偏窄的是第三種做法。 只用每人 3 次的平均去算差值、然後不做任何修正，跨距是 74.22 mmHg，窄了 9.62%，少掉 7.90 mmHg。這才是 Bland 與 Altman 1999 那篇的第 5.2 節在警告的東西：平均會把每個方法自己的測量誤差壓小，於是你報出來的界限描述的是「三次平均對三次平均」，而臨床上每個病人只會被量一次。

為什麼界限沒有跟著垮

因為 J 與 S 是同時讀同一位受試者的。那一瞬間血壓真正的高低，兩個方法共享，於是它在同一次的差值裡被減掉了。差值的變異因此已經幾乎不含「這個人這一刻血壓在哪裡」的成分，重複三次也不會憑空多出來——把它們當成獨立觀測，估到的 SD 自然與正確答案接近。

正確的做法是把差值的變異拆開再組回去。以 $m$ 表示每人的重複次數、$\bar d$ 表示每人的平均差值、$s_w^2$ 表示各方法自己的組內變異：

代進這份資料：$m$ 是 3，修正係數 0.667；每人平均差值的標準差 18.93 mmHg；觀察者 J 自己重複三次的組內標準差 6.12 mmHg，機器 S 是 9.12 mmHg（各自 170 個自由度）。組回去得到 20.95 mmHg，就是上表 C 那一橫列。

順帶注意兩個組內標準差的大小關係：機器自己重複三次的變異比人工還大（9.12 對 6.12 mmHg）。一台號稱能消除觀察者誤差的裝置，在自身再現性上並沒有比人穩。

比例性偏誤：檢查過，未偵測到

一組固定的界限有一個隱含前提：差值的大小與測量值的大小無關。血壓高的人差得比較多，就不成立——這時候一條界限套整個量程，會在低端過寬、在高端過窄。

檢查方式是把差值對平均值做線性迴歸，看斜率。這份資料上，以每人平均計算的斜率是 +0.0356 mmHg/mmHg（95% CI −0.1028 到 +0.1739，p = 0.610）；換成全部 255 個單次配對是 +0.0514（p = 0.226）。點估計的意思是測量值每上升 10 mmHg，差值變動 +0.36 mmHg——與 10 mmHg 這個量級相比可以忽略。

所以本頁的結論是：在這份資料裡未偵測到比例性偏誤，界限就用一組固定的數字報。

右邊那個面板示範另一種常見做法：兩邊取對數再做同樣的事。這在差值隨測量值成比例放大時是標準處理，好處是界限會以比值而不是絕對量回來。這份資料上取完對數斜率是 −0.0640（p = 0.338），一樣未偵測到趨勢；對數尺度的結果換算回來是：機器平均讀高 12.6%，界限為比值 0.851 到 1.491，也就是 -14.9% 到 +49.1%。

動手跑一次

data(sbp, package = "MethComp")

# 資料是長格式：meth（J/R/S）、item（受試者）、repl（第幾次）、y（mmHg）。
# 先確認每個 (meth, item, repl) 只有一列，否則下面的 reshape 會靜默地取到
# 第一筆而不報錯——這一行是承重的，不是裝飾。
stopifnot(!anyDuplicated(sbp[c("meth", "item", "repl")]))

w <- reshape(sbp[sbp$meth %in% c("J", "S"), ],
             direction = "wide", idvar = c("item", "repl"),
             timevar = "meth")
w$diff <- w$y.S - w$y.J          # 差值定義：機器減人工
w$avg  <- (w$y.S + w$y.J) / 2

# 每人的平均（85 位受試者）
m <- aggregate(cbind(diff, avg, y.J, y.S) ~ item, data = w, FUN = mean)

# 1. 誘導：相關係數。高，而且與一致性無關。
cor.test(m$y.J, m$y.S)

# 2. 一致性界限，Bland & Altman (1999) 的重複測量修正。
#    直接對 255 個配對算 sd(w$diff) 會把 85 個人當成 255 個獨立觀測：
#    界限本身只窄百分之幾，但界限的信賴區間會窄三分之一。
mm <- 3
within_var <- function(v) {
  # 每個方法自己的組內變異：各受試者的變異數，用自由度加權平均
  a <- aggregate(v ~ item, data = w, FUN = function(x) c(var(x), length(x) - 1))
  sum(a$v[, 1] * a$v[, 2]) / sum(a$v[, 2])
}
n      <- nrow(m)
df_w   <- n * (mm - 1)
k      <- 1 - 1 / mm
vJ     <- within_var(w$y.J)
vS     <- within_var(w$y.S)
v_dbar <- var(m$diff)
s2     <- v_dbar + k * (vJ + vS)
bias   <- mean(m$diff)
z      <- qnorm(0.975)
loa    <- bias + c(-1, 1) * z * sqrt(s2)

# 3. 界限的信賴區間。s2 是三個獨立成分加起來的，它的變異數也就是三份
#    Var(s2) = 2 s2^2 / df 加起來；界限的變異數再由
#    Var(LoA) = Var(bias) + z^2 * Var(s) 組出來，其中 Var(s) = Var(s2) / (4 s2)。
#    課本上那條 sqrt(s2 * (1/n + z^2 / (2(n-1)))) 只在「所有觀測獨立」時成立，
#    這裡不成立：套下去界限的區間會比上表寬大約一個 mmHg。分母是受試者數，
#    不是配對數。
var_s2 <- 2 * v_dbar^2 / (n - 1) + k^2 * 2 * vJ^2 / df_w + k^2 * 2 * vS^2 / df_w
se     <- sqrt(v_dbar / n + z^2 * var_s2 / (4 * s2))
tq     <- qt(0.975, n - 1)
rbind(lower = loa[1] + c(-1, 1) * tq * se,
      upper = loa[2] + c(-1, 1) * tq * se)

# 4. 比例性偏誤：差值對平均值迴歸。橫軸必須是兩者的平均，
#    對參考方法作圖會憑空生出一條向下的斜率。
summary(lm(diff ~ avg, data = m))

# 5. 對數尺度版：界限回來時是比值，不是 mmHg。跟第 4 步一樣跑在每人平均上——
#    跑在 255 個單次配對會得到另一個斜率，那不是上表報的那個數字。
lm(log(y.S) - log(y.J) ~ I((log(y.S) + log(y.J)) / 2), data = m)

# 匯出給下面的 Python 用
write.csv(sbp, "sbp.csv", row.names = FALSE)

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats

# MethComp 沒有登記在 Rdatasets 上，get_rdataset("sbp", "MethComp") 會回 404。
# 上面那段 R 已經把它寫出來了。
sbp = pd.read_csv("sbp.csv")
assert not sbp.duplicated(["meth", "item", "repl"]).any()

w = (sbp[sbp.meth.isin(["J", "S"])]
     .pivot(index=["item", "repl"], columns="meth", values="y")
     .reset_index())
w["diff"] = w["S"] - w["J"]
w["avg"] = (w["S"] + w["J"]) / 2
m = w.groupby("item")[["diff", "avg", "J", "S"]].mean()

print(stats.pearsonr(m["J"], m["S"]))

mm = 3
def within_var(col):
    g = w.groupby("item")[col]
    return (g.var(ddof=1) * (g.count() - 1)).sum() / (g.count() - 1).sum()

n = len(m)
df_w = n * (mm - 1)
k = 1 - 1 / mm
vJ, vS = within_var("J"), within_var("S")
v_dbar = m["diff"].var(ddof=1)
s2 = v_dbar + k * (vJ + vS)
bias = m["diff"].mean()
z = stats.norm.ppf(0.975)
loa = bias + np.array([-1, 1]) * z * np.sqrt(s2)

var_s2 = 2 * v_dbar**2 / (n - 1) + k**2 * 2 * vJ**2 / df_w + k**2 * 2 * vS**2 / df_w
se = np.sqrt(v_dbar / n + z**2 * var_s2 / (4 * s2))
tq = stats.t.ppf(0.975, n - 1)
print(loa, loa[:, None] + np.array([-1, 1]) * tq * se)

print(smf.ols("diff ~ avg", data=m).fit().summary())

lg = m.assign(ld=np.log(m["S"]) - np.log(m["J"]),
              la=(np.log(m["S"]) + np.log(m["J"])) / 2)
print(smf.ols("ld ~ la", data=lg).fit().params)

怎麼讀報表

一篇方法比較的論文，Results 通常給你一張 Bland-Altman 圖加一段文字。六個位置值得逐一確認：

1. 差值的方向寫清楚了嗎。 「A 減 B」與「B 減 A」的偏差正負號相反。圖的縱軸標籤應該寫得出是誰減誰；只寫 “difference” 的圖沒辦法被引用。
2. 界限有沒有附信賴區間。 界限是估計值。本頁的上界是 +56.7，它的 95% CI 一路到 +63.3——若作者要主張「界限落在可接受範圍」，要看的是區間的外緣，不是點估計。
3. 重複測量怎麼處理的。 有重複而 n 報成配對數（本頁的話是 255 而不是 85），那篇的界限信賴區間偏窄；只用每人平均而不修正的，界限本身就偏窄。Methods 沒交代的，兩種都有可能。
4. 比例性偏誤檢查了嗎，怎麼檢查的。 應該看得到「差值對兩者平均」的迴歸斜率。若是對參考方法迴歸的，那條斜率不能信。
5. 有沒有給臨床上可接受的差，而且是事先給的。 沒有這個數字，「界限可接受」就只是作者的感想。最好的版本是引用外部標準，像本頁的 ISO 81060-2 criterion 1。
6. 樣本涵蓋的量程夠不夠。 界限只能外推到有資料的範圍。血壓集中在正常區間的一群人，說不出這台機器在 140 mmHg 以上準不準。

常見誤用

這一頁與其他頁的關係

• 同一份資料的信度那一面——見組內相關係數（ICC）。那一頁用的是本頁刻意避開的 J 對 R，兩頁合起來就是「一致性」與「信度」是同一批數字上的兩個不同問題。
• 類別變項的一致性——見Kappa 與加權 Kappa。結果是類別而不是連續量時，差值沒有意義，換一套工具。
• 測量誤差的下游後果——見測量誤差。本頁量到的差值分散度，就是那一頁在講的東西進到迴歸之前的樣子。
• 差值對平均值的那條迴歸——見線性迴歸。比例性偏誤的檢查只是把它用在一個特別的 xy 上。
• 重複測量的相依結構——見叢集資料。本頁的修正公式是那一頁變異數分解的一個特例。
• 切點與決策翻轉——見切點的選擇。本頁第二條翻譯路徑用的就是那一頁的框架。
• 臨床脈絡——見設計章 診斷準確度研究。

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B9-03-bland-altman.R

回報這一頁的問題