R包npmc的非参数多重比较

单因素方差分析（单因素ANOVA）后接多重比较（Tukey HSD检验），是比较多组数据两两差异的常规做法。

不过由于单因素ANOVA的前提假设条件比较严格，要求数据必须满足正态性、方差齐性等，而很多情况下我们的数据并不符合ANOVA的条件，因此多重比较也无法直接使用。此时除了考虑转换数据外，另一种思路就是使用非参数的方法。

对于非参数的方法，此前提到过可以首先执行Kruskal-Wallis检验代替单因素ANOVA比较多组数据的整体差异，然后再通过Wilcoxon检验代替Tukey HSD继续确定两组间差异水平，并通过p值校正控制I类错误率。

另一种常见做法则是使用非参数的多重比较。《R语言实战 第一版》（154页）提到了一个R包，npmc，提供了一种非参数多组比较程序npmc()，可在控制I类错误的前提下，进行所有组之间的成对比较。

注：《R语言实战 第二版》同位置处就更改为wilcoxon检验+p值校正的方法了，不是非参数多组比较程序本身方法的问题，而是npmc包在R3.0以上的版本不再支持。因此若期望使用该程序包，还需安装一个较早版本的R（低于3.1版本）。（这也是为什么我此前并不知道npmc包方法的原因，前文也吐槽过这一点）

接下来，展示使用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')

友情链接