R包npmc的非引數多重比較

2025-04-29 07:35 小白魚的生統筆記

R包npmc的非引數多重比較
不過由於單因素ANOVA的前提假設條件比較嚴格，要求資料必須滿足正態性、方差齊性等，而很多情況下我們的資料並不符合ANOVA的條件，因此多重比較也無法直接使用。此時除了考慮轉換資料外，另一種思路就是使用非引數的方法。
另一種常見做法則是使用非引數的多重比較。《R語言實戰第一版》（154頁）提到了一個R包，npmc，提供了一種非引數多組比較程式npmc()，可在控制I類錯誤的前提下，進行所有組之間的成對比較。
接下來，展示使用npmc包的方法，完成非引數的多重比較過程。
示例資料集
npmc的內建資料集brain，來自某項對大腦疾病的研究，左半球（LH）或右半球（RH）病變對反應時間的影響。
library(npmc)
#資料集，詳情 ?brain

data(brain)

head(brain)

受試者包括30名對照（c）、30名LH腦損傷的患者（l）和30名 RH腦損傷的患者（r），是一個平衡的單因素研究案例。
接下來期望比較正常人群、LH或RH腦損傷的患者，在反應時間上是否存在顯著不同。
可能首先會想到使用單因素ANOVA來做。
但由於原始資料不滿足方差齊性，無法透過單因素ANOVA解決，所以考慮使用非引數的方法。
#Bartlett 檢驗顯示 p 值大於 0.05，即拒絕方差齊性假設

bartlett.test(var~class, data = brain)

npmc的非引數多重比較
在進行非引數多重比較之前，可首先使用Kruskal-Wallis檢驗檢視整體差異。
#Kruskal-Wallis 檢驗，整體差異

fit <- kruskal.test(var~class, data = brain)

fit$p.value

結果顯示3組資料在整體上存在顯著差異。
接下來，繼續使用npmc包中的方法進行兩組間的比較。
npmc()函式中提供了兩種非引數多重比較方法，Behrens-Fisher和Steel，均可用於成對和多對一的比較情況，並同時計算相對效應的置信區間（1-alpha）。
npmc()函式接受的輸入為一個兩列的資料框，其中一列名為class（分組變數），另一列名為var（因變數）。因此在輸入資料前，更改變數名稱保證能夠被npmc()識別。
npmc包的功能實現需要mvtnorm包計算檢驗統計量的p值，如果mvtnorm包不可用，則結果將包含NA的p值。mvtnorm包中的函式使用隨機值進行積分計算，因此npmc()關於p值和置信區間的結果因呼叫而異，並且只能被視為近似解。
#npmc 的非引數多重比較

npmc_test <- npmc(brain, df = 2, alpha = 0.05)
#概要

summary(npmc_test, type = 'both', short = FALSE)

#或者

npmc_test$test

對於結果輸出概要的解釋：
BF/Steel，分別為Behrens-Fisher和Steel的檢驗結果；
cmp，比較的組名；
gn，兩組樣本量之和；
effect，相對效應估計值；
variance，方差估計值；
std，標準偏差；
statistic，檢驗統計量；
p-value 1s，單側檢驗p值，與雙側檢驗相比它只考慮一端情況，若a-b比較組中a和b存在顯著差異則代表a顯著小於b；
p-value 2s，雙側檢驗p值；
zero，TRUE代表0方差或近似0方差，否則為FALSE。
一般情況下，關注雙側檢驗的結果就可以了。
對於本示例，儘管Kruskal-Wallis檢驗顯示資料整體存在顯著差異，但p值比較接近0.05顯著性閾值了，暗示兩組間的差異可能不明顯。
非引數多重比較結果中，Behrens-Fisher顯示l組（LH腦損傷的患者）和r組（RH腦損傷的患者）的反應時間存在區別，但不明顯；Steel則顯示無差異。
最後畫個箱線圖描述資料分佈特徵。
#ggplot2 箱線圖

library(ggplot2)
p <- ggplot(data = brain, aes(x = class, y = var, fill = class)) +

geom_boxplot(outlier.size = 1) +

theme(panel.grid = element_blank(), panel.background = element_rect(fill = 'transparent', color = 'black'),

legend.title = element_blank(), legend.key = element_blank(), plot.title = element_text(hjust = 0.5)) +

labs(x = '', y = '', title = 'Behrens-Fisher\n')
p +

annotate('segment', x = 1, xend = 2, y = 50, yend = 50) +

annotate('segment', x = 1, xend = 3, y = 55, yend = 55) +

annotate('segment', x = 2, xend = 3, y = 60, yend = 60) +

annotate('text', x = 1.5, y = 52, label = 'p = 0.147') +

annotate('text', x = 2, y = 57, label = 'p = 0.754') +

annotate('text', x = 2.5, y = 62, label = 'p = 0.044')

友情連結

相關文章