無母數檢定：Wilcoxon 與 Kruskal-Wallis

這個方法在解決什麼問題

「資料看起來不常態，所以我改跑 Wilcoxon」——這是無母數檢定最常見的登場方式， 而它幾乎總是伴隨一個沒有被說出口的誤解：以為 Wilcoxon 是在比中位數。

它不是。Wilcoxon rank-sum 檢定（又叫 Mann-Whitney U 檢定）比的是隨機優勢 （stochastic ordering）：隨機抽一個 A 組的人、隨機抽一個 B 組的人，前者比後者大的機率是不是二分之一。 寫成式子就是虛無假設 $P(X > Y) = 1/2$。中位數只是分布的一個摘要，而這個機率吃的是整個分布。

兩者大部分時候會給出一致的結論，但不是同一件事——而它們分開的時候， 報表上就會出現一組看起來自相矛盾的數字：兩組中位數印在 Table 1 上一模一樣， 底下那個 p 值卻是顯著的。本頁的第二節就是把那個情境完整做出來。

本頁用兩份資料：MASS::birthwt（189 位產婦，結果變項是新生兒出生體重，公克） 與 datasets::sleep（Student 1908 年那份安眠藥資料）， 與 t 檢定與變異數分析那一頁用的是同一組變項定義， 所以兩頁的數字可以直接對照。

虛無假設是隨機優勢，不是中位數

下面這組資料是刻意構造出來的，用來示範一件真的會發生的事。 情境是同一種手術在兩家醫院的住院天數（單位是天）：

• 甲院沒有固定的出院流程，每個病人依狀況決定，於是有一條很長的下尾（術後隔天就走的人不少），上方則封在 9 天。
• 乙院有一套以第七天出院為目標的流程，所以幾乎沒有人在第 6 天以前離開，剩下的變異全部往上跑——併發症把少數人留到三、四週。

甲院（30 人）：1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9

乙院（30 人）：6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 9, 10, 11, 12, 13, 14, 17, 18, 20, 22, 25

兩組中位數的差是 0 天——完全相同。 而 Wilcoxon rank-sum 檢定給出 W = 629，p = 0.007。

中間那格是重點。乙院的累積分布函數從頭到尾都不高於甲院—— 也就是說，不論你把門檻設在幾天，乙院超過那個門檻的比例都不會比甲院低。 這正是「乙院系統性地比較久」的嚴格意思，而它與「兩組中位數相同」完全相容： 兩條階梯都在第 7 天那一階跨過 0.5。

順帶一提，同一組資料的平均數差是 3.87 天，Welch t 檢定 p < 0.001。 平均數看得見的東西，中位數看不見——右尾那幾個住到三週的病人，正是臨床上最該被看見的人。

把隨機優勢變成一個能報的數字

$P(X > Y)$ 不必停在抽象的層次，它可以直接估出來：把所有跨組配對拿出來數一遍， 一方較大就記一分，兩者相同就各記半分，除以配對總數。上圖的 C 格畫的就是這件事—— 900 個格子，每一格是一對病人。

構造例的估計值是 0.699，也就是說隨便抓一組甲乙兩院的病人， 乙院那位住得比較久的機率大約是 70%。

真實資料也一樣可以報。birthwt 裡未抽菸與抽菸兩組的出生體重， 這個機率是 0.617——隨機抓一對母嬰，未抽菸那組的嬰兒比較重的機率大約 62%。

它的好處是可以直接講給人聽：「隨機挑一位接受新療法的病人，他的恢復比隨機挑一位接受標準治療的病人好， 機率是七成」，比 p 值有訊息量得多，而且不需要任何分布假設。

只報 p 值不夠：Hodges-Lehmann 位移估計

Wilcoxon 的輸出預設只有一個 p 值，於是很多論文的 Results 就寫成 「中位數 A vs B（p = …）」。這個寫法把兩個不同的東西併在一起： 中位數是各組的摘要，p 值是組間的檢定，中間缺了「差多少、精確度多少」的效果量。

補上那一格的是 Hodges-Lehmann 位移估計：所有跨組配對差值的中位數。 R 的 wilcox.test(conf.int = TRUE) 會一併給出它與信賴區間，單位與原始資料相同。

以 birthwt 的抽菸與否為例，同一個比較有三個不同的數字：

三個數字都在三百公克上下，但它們不是同一個量。中間那列的定義值得再讀一次： Hodges-Lehmann 估的是「配對差值的中位數」，不是「中位數的差」—— R 的說明文件把後者明白標成常見誤解。這兩者在對稱分布下會很接近，在偏態分布下不會。

實務上的寫法：報「Hodges-Lehmann 位移估計 306.2 公克 （95% CI 85 到 512，Wilcoxon rank-sum p = 0.007）」， 必要時再加上 $P(X > Y)$ 的估計 0.617。 順帶檢查：直接算 8510 個配對差值的中位數是 307 公克， 與 wilcox.test() 回報的 306.2 公克非常接近，但不是同一個數。 差別來自算法：wilcox.test() 是把檢定反轉、用 uniroot 解出來的， 不是直接取那 8510 個差值的中位數。

無母數還是取對數？兩條路估的對象不同

右偏的變項（住院天數、費用、生物標記濃度）還有第二條路：取對數之後跑 t 檢定。 兩條路都合理，但它們的產出不是同一個量，而論文必須說清楚走的是哪一條。

同一份 birthwt 資料，兩條路的結果是：

第二列的意思是：未抽菸組的幾何平均出生體重（2950.8 公克） 是抽菸組（2679.8 公克）的 1.101 倍， 換算成百分比是高出 10.1%（95% CI 1.4% 到 19.6%）， p = 0.022。

差別不只是換算單位而已：

• 加法的位移假設兩組相差一個固定的量。出生體重適合這樣看——臨床門檻（低出生體重 2500 公克）是加法的，「少 300 公克」對每一個嬰兒的意義大致相同。
• 乘法的比值假設兩組相差一個固定的倍數。跨越好幾個數量級的量（病毒量、C 反應蛋白、住院費用）適合這樣看——「高 10%」在低值與高值處都成立，而「高 300 單位」在低值處可能是三倍、在高值處微不足道。
• 取對數的 t 檢定估的是幾何平均的比值，不是算術平均的比值。把對數尺度上的平均數反轉回來得到的是幾何平均，它一定小於或等於算術平均；要主張算術平均要另外處理。
• 取對數要求所有觀測值為正。有零值時常見的 log(x + 1) 會讓那個常數影響結果，而且比值失去乾淨的解釋。

相同值與精確 p 值

秩檢定的精確 p 值來自「把 N 個觀測值重新分派到兩組」的所有排列。 資料裡出現相同值（ties）時，這些觀測值只能給相同的平均秩，排列的參考分布因此改變。

birthwt 的出生體重記錄到公克，仍然有 131 個相異值分給 189 筆資料， 其中 97 筆與別人共用同一個值，最大的一群有 5 筆。 臨床資料幾乎不可能沒有 ties。

多數教科書寫的是「有 ties 就算不出精確 p 值，軟體會退回常態近似」。 這句話對現在的 R 已經不成立（本頁驗證於 R 4.6.0）： wilcox.test() 內建了 Streitberg-Röhmel 移位演算法， 在有 ties 時改用「條件於觀測到的秩」的排列分布，照樣給精確推論。

真正切換預設的是樣本數，不是 ties：兩組都少於 50 筆時走精確法，否則走常態近似。 所以本頁的構造例（每組 30 筆）拿到的是精確 p 值， 而 birthwt 的抽菸比較（115 對 74）拿到的是常態近似。

把 birthwt 裡出生體重在全體中唯一的那些列抽出來（92 筆，完全沒有 ties）， 再把同一批資料四捨五入到最接近的 250 公克—— 也就是病歷回顧實際上能拿到的精度——可以看到 ties 的代價：

同一批病人、同一個比較，只因為量得比較粗，檢定統計量與 p 值就都動了。 測量精度是分析的一部分，不是資料收集階段結束後就不再影響結果的事。

幾個實務要點：

• 報軟體與版本。 同一份資料在舊版 R、SPSS、SAS 上可能拿到常態近似而不是精確 p 值，兩者不會完全相同。
• 連續性校正只影響常態近似那條路。 本頁的抽菸比較有校正時 p = 0.00677，沒有校正時 p = 0.00674；差距通常很小，但 p 值貼著門檻時就不是小事。
• 不要為了消除 ties 去加隨機擾動。 那等於把答案交給亂數種子。
• ties 多到某個程度時，該換的是模型不是檢定：結果變項本來就是有序類別（疼痛分級、NYHA 分級）時，有序邏輯迴歸能同時給出效果量與共變項校正。

配對版：Wilcoxon signed-rank

同一個人前後測、同一個人左右眼、配對設計——有一對一的連結就要用配對版本， 理由與 配對 t 檢定完全相同：配對把人與人之間的差異整個消掉。

Wilcoxon signed-rank 檢定的做法是：算出每個人的差值，依絕對值大小排序給秩， 再看正差值的秩和有沒有偏離。用 Student 那份安眠藥資料（10 位受試者各試兩種藥）， 每個人的差值是 1.2, 2.4, 1.3, 1.3, 0, 1, 1.8, 0.8, 4.6, 1.4。

注意其中有 1 個差值為零的受試者。這在配對資料裡很常見， 而不同軟體的處理方式不一樣：

兩種零差值處理法在這份資料上給出相同的 p 值，但點估計與區間都不一樣， 而且沒有任何輸出會提醒你剛剛做了一個選擇。有零差值時，把軟體與版本寫進 Methods。

還有一個容易被忽略的細節：上表兩個無母數區間的實際涵蓋水準是 95.6% 與 96.5%，不是剛好 95%。 秩統計量是離散的，能達到的水準只有有限幾個，R 會挑最接近的並如實標出來。

三組以上：Kruskal-Wallis 與事後兩兩比較

Kruskal-Wallis 檢定是 ANOVA 的秩版本：把所有觀測值一起排序，再問各組的平均秩是不是都一樣。 它同樣是一個綜合檢定（omnibus test），顯著只代表「不全相同」，不說是誰跟誰不同。

用 birthwt 依產婦族群分成三組：

綜合檢定的結果是 $\chi^2$ = 8.52，df = 2，p = 0.014。

接下來的兩兩比較有一個常見的做法是錯的：跑三次 Wilcoxon rank-sum。 問題在於每一次都只拿其中兩組重新排秩，第三組的資訊被丟掉了， 而且各次的參考分布互不相容。正確的做法（Dunn 檢定）是只排一次秩， 用全體樣本算出的變異數去比較各組的平均秩，並且共用同一套 ties 校正。

這張表是本節的重點。綜合檢定 p = 0.014 是顯著的， 但三個兩兩比較經 Holm 調整之後沒有一個達到統計顯著。 這不是矛盾，也不代表各組之間沒有差異——只代表這份資料的證據強度不足以在控制家族錯誤率的前提下， 指認出是哪兩組不同。要寫的是「未達統計顯著」，不是「沒有差異」。

幾個實務要點：

• 事後比較的方法與調整法要事先寫進計畫書。 「所有兩兩比較」與「只比對照組」需要的調整不同，型一錯誤累積的原理見 t 檢定與變異數分析。
• Holm 幾乎總是優於 Bonferroni。 兩者控制的都是家族錯誤率，Holm 是逐步做法，檢定力嚴格較高而假設不變。
• Kruskal-Wallis 也不是在比中位數。 要把它讀成「位移」需要各組分布形狀相似這個額外假設，和兩組的情況一樣。
• 有共變項要校正時，秩檢定就不夠用了。 秩檢定沒有校正機制；改用有序邏輯迴歸或對轉換後的變項跑線性迴歸。

動手跑一次

library(MASS)
data(birthwt, package = "MASS")

birthwt$smoke_f <- factor(birthwt$smoke, levels = c(0, 1),
                          labels = c("Non-smoker", "Smoker"))
birthwt$race_f  <- factor(birthwt$race, levels = 1:3,
                          labels = c("White", "Black", "Other"))

# 兩組獨立樣本：一定要加 conf.int = TRUE，才拿得到 Hodges-Lehmann 估計
wilcox.test(bwt ~ smoke_f, data = birthwt, conf.int = TRUE)

# P(X > Y)：直接數配對，同值各算半分
ns <- birthwt$bwt[birthwt$smoke_f == "Non-smoker"]
sm <- birthwt$bwt[birthwt$smoke_f == "Smoker"]
mean(outer(ns, sm, ">")) + 0.5 * mean(outer(ns, sm, "=="))

# 另一條路：取對數後 t 檢定，得到的是幾何平均的比值
fit <- t.test(log(ns), log(sm))
unname(exp(c(diff(rev(fit$estimate)), fit$conf.int)))

# 配對版。⚠️ R 4.6 起公式介面不接受 paired = TRUE（會直接報錯），
# 要嘛給兩個向量，要嘛用 Pair() 的寫法。
data(sleep)
d1 <- sleep$extra[sleep$group == 1]
d2 <- sleep$extra[sleep$group == 2]
wilcox.test(d2, d1, paired = TRUE, conf.int = TRUE)
wilcox.test(Pair(d2, d1) ~ 1, conf.int = TRUE)   # 等價寫法

# 三組以上
kruskal.test(bwt ~ race_f, data = birthwt)

# Dunn 事後比較：只排一次秩，用全體的變異數，最後做 Holm 調整
r <- rank(birthwt$bwt); N <- length(r); tie <- table(birthwt$bwt)
s2 <- N * (N + 1) / 12 - sum(tie^3 - tie) / (12 * (N - 1))
nn <- tapply(r, birthwt$race_f, length); rb <- tapply(r, birthwt$race_f, mean)
cb <- combn(levels(birthwt$race_f), 2)
z  <- apply(cb, 2, function(k) (rb[[k[1]]] - rb[[k[2]]]) /
                               sqrt(s2 * (1 / nn[[k[1]]] + 1 / nn[[k[2]]])))
data.frame(pair = apply(cb, 2, paste, collapse = " vs "),
           z = round(z, 2),
           p_holm = p.adjust(2 * pnorm(-abs(z)), method = "holm"))

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

bw = sm.datasets.get_rdataset("birthwt", "MASS").data
ns = bw.loc[bw["smoke"] == 0, "bwt"].to_numpy()
smk = bw.loc[bw["smoke"] == 1, "bwt"].to_numpy()

# 檢定本身有：method="exact" / "asymptotic" 對應 R 的 exact 參數
stats.mannwhitneyu(ns, smk, alternative="two-sided")

# ⚠️ scipy 沒有 Hodges-Lehmann 估計，也沒有它的信賴區間。
# 點估計可以自己算（就是配對差值的中位數）；區間要自己 bootstrap，
# 那與 R 用檢定反轉得到的區間不是同一個東西，不要混報。
diffs = np.subtract.outer(ns, smk).ravel()
hodges_lehmann = np.median(diffs)

# P(X > Y) 同樣自己數
p_superior = (np.mean(np.subtract.outer(ns, smk) > 0)
              + 0.5 * np.mean(np.subtract.outer(ns, smk) == 0))

sleep = sm.datasets.get_rdataset("sleep").data
d1 = sleep.loc[sleep["group"] == 1, "extra"].to_numpy()
d2 = sleep.loc[sleep["group"] == 2, "extra"].to_numpy()
stats.wilcoxon(d2, d1)          # zero_method 預設與 R 不同，要顯式指定

groups = [g["bwt"].to_numpy() for _, g in bw.groupby("race")]
stats.kruskal(*groups)
# ⚠️ Dunn 事後比較不在 scipy 裡（scikit-posthocs 有）；
# 要留在標準函式庫裡就照上面 R 那段的公式自己寫。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B1-04-nonparametric.R

延伸觀看

▶醫學統計 EP10 t 檢定與非參數法

▶醫學統計 EP11 ANOVA

回報這一頁的問題