【原】R笔记：多个独立样本的非参数检验及其多重比较

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：多个独立样本的非参数检验及其多重比较。

方差分析是假设数据满足正态分布和方差齐性的前提下，计算总体相等但在当前样本量下由于偶然原因导致均值出现差异的概率，或者说基于当前样本量，结局均值的差异有多大的可能性（概率）是由于偶然原因造成的。若不满足前提条件则宜采用更为稳健的非参数方法，比如常用的基于秩次的假设检验（秩和检验）。秩和检验时，我们的H0假设不再是均值相等，而是中位数代表的分布是否相同。

类似于方差分析中的独立样本的方差分析和区组设计的方差分析，多个样本的非参数检验分为多个独立样本的非参数检验和多个相关样本的非参数检验。此次笔记是关于多个独立样本的秩和检验及其事后两两比较方法，示例1是结局变量为连续变量的例子，示例2则是解释/分组变量为分类变量，结局/指标变量为有序分类变量。

此次笔记涉及到分析方法是单因素方分析、Kruskal-Wallis检验及事后两两比较、Ridit分析等，涉及命令有shapiro.test, leveneTest, oneway.test, kruskal.test, posthoc.kruskal.nemenyi.test, posthoc.kruskal.dunn.test, ridit等。

示例1：结局变量为连续变量。

示例1

正常、单纯性肥胖及皮质醇增多症三组人的血浆皮质醇含量(nmol/L)的差异有无统计学意义？

多个独立样本，连续型资料，可以考虑的是方差分析。方差分析需要满足正态分布和方差齐性。library(readxl)

inpt <- read_excel("D:/Temp/NPtdata.xlsx")  #从excel中导入数据

by(inpt$cortisol,inpt$G,shapiro.test)  #正态性检验

library(car) 

leveneTest(inpt$cortisol,inpt$G)  #方差齐性检验 

本例可以采用单因素的方差分析，可参照我们在《R笔记：单因素方差分析》一文中介绍的aov {stats}，lm {stats}。本例另外一个函数oneway.test{stats}，在正态分布下，不论方差相同与否均可使用该函数：oneway.test(formula, data, subset, na.action, var.equal = FALSE)oneway.test(cortisol~G,data=inpt,var.equal = TRUE) #单因素方差分析

当然，本次笔记的重点是非参数检验的操作演示。当数据不满足正态分布或方差齐性时，需要采用更为稳健的分析方法，比如非参数法。当然，满足方差分析的条件也是可以使用非参数方法的，只是检验效能低一点而已。本例采用Kruskal-Wallis H检验。Usage:

kruskal.test(x, g, ...)

kruskal.test(formula, data, subset, na.action, ...)

本例命令：

kruskal.test(inpt$cortisol,inpt$G)

或

kruskal.test(cortisol~G,data=inpt)

当总体检验有统计学意义后，接下来我们会想知道哪两个组间会存在差异，就涉及到两两比较（或者叫多重比较）。推荐程序包PMCMR(The Pairwise Multiple Comparison of Mean Ranks Package)，为非默认包，使用前需要先下载[install.packages("PMCMR")]。该程序包提供了大量的秩和检验（包括Kruskal Wallis、Friedman）后的多重比较方法，如下。我们仅以常见的Nemenyi、Dunn来进行演示。

①Nemenyi法：posthoc.kruskal.nemenyi.test {PMCMR}

Usage:

posthoc.kruskal.nemenyi.test( x, g, dist =c("Tukey", "Chisquare"), ...)

posthoc.kruskal.nemenyi.test(formula, data, subset,na.action, dist =c("Tukey", "Chisquare"), ...)

library(PMCMR)

posthoc.kruskal.nemenyi.test(cortisol~G,data=inpt)结果显示单纯肥胖组和正常组的血浆皮质醇无统计学差异，但皮质醇增多症组与单纯肥胖组(P=0.002)、皮质醇增多症组与正常组（P<0.001）均有统计学差异。但结果也显示，数据存在“结”（重复数据），需要进行校正。

在Nemenyi法中，Only for method = "Chisq" a tie correction is employed.

posthoc.kruskal.nemenyi.test(inpt$cortisol,inpt$G,dist="Chisquare")经校正后，P值有所变化，但最终结论未变：皮质醇增多症组与单纯肥胖组(P=0.003)、皮质醇增多症组与正常组（P<0.001）均有统计学差异。最后的提示也表明采用了卡方对“结”进行了校正。

②Dunn法：posthoc.kruskal.dunn.test {PMCMR}Usage:

posthoc.kruskal.dunn.test( x, g, p.adjust.method =p.adjust.methods, ...)

posthoc.kruskal.dunn.test(formula, data, subset,na.action, p.adjust.method = p.adjust.methods, ...)library(PMCMR)

posthoc.kruskal.dunn.test(inpt$cortisol,inpt$G)

posthoc.kruskal.dunn.test(cortisol~G,data=inpt, p.adjust.method ="bonferroni") #数据存在“结”，采用bonferroni进行校正。P值校正方法有"holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none"数据存在“结”时，dunn法会直接给与校正，默认给与了holm校正。结果解读可参照Nemenyi法，也不再累述。

PMCMR加强版为PMCMRplus，提供了更为丰富的比较方法，下载安装[install.packages("PMCMRplus")]，不再累述。示例2：解释/分组变量为分类变量，结局/指标变量为有序分类变量。

示例2

244例疥疮患者接受2.5%和1%的二硫化硒洗剂、10%的硫磺软膏治疗，其中有效药物硫磺软膏外用作为对照，考察二硫化硒治疗疥疮的疗效。

岭南皮肤性病科杂志,2005,12(1):57-58.

在很多杂志期刊中，此类数据常误用卡方来进行分析。这种分组变量为有序/无序分类变量、指标变量为有序分类变量的资料，应该使用非参数检验法。三组独立样本，宜采用前面刚刚讲过的Kruskal-Wallis检验。文中使用的并不是秩和检验法，而是Ridit法，我们也会在此演示一下。library(readxl)

inptc <- read_excel("D:/Temp/NPtdata.xlsx",2)  #从excel中导入第2个sheet的数据，sheet2中的数据是原始的数据格式。笔者为找到R中Kruskal-Wallis检验的数据频数加权方法，需要分析的数据需要是原始的数据格式（指按照每条记录占一行的原始格式，而不是概括成频数的格式，我们在《R笔记：两独立样本的差异显著性检验》已有过说明）

kruskal.test(inptc$eft,inptc$group) #或kruskal.test(eft~group,data=inptc)结果显示：3个治疗组的疗效无统计学差异（chi2=2.780，P=0.249）。如果结果有统计学意义，可以再进行两两比较，本例总体比较无统计学意义，两两比较也就失去了意义，本例仅作演示。

library(PMCMR)

posthoc.kruskal.dunn.test(eft~group,data=inptc)

秩和检验时最为常用的非参数检验方法，但不是唯一，比如这种结局变量是等级资料的数据，我们也可以采用Ridit分析（relative to an identified distribution unit）。Ridit分析是将有序分类资料（等级资料）转换成成连连续数据（即Ridit值），然后可按照正态分布的理论对Ridit值进行比较。在转换时需要先确定一个参照组或者叫标准组，确定参照组后可获得参照组各等级的Ridit值及平均Ridit值（平均Ridit恒=0.5），其他各组的平均Ridit需要根据参照组各等级的R值来计算，值在0-1之间。参照组的确定一种方法是所有组的总和，另外一种是选观察人数较多、数据比较稳定的组。

ridit {Ridit}：An extension of the Kruskal-Wallis TestUsage:

ridit(x, g, ref = NULL)

x: a numeric vector of data values or a matrix of crosstab data.

g: a vector giving group of data or when x is a crosstab, number 1 or 2 when group is in the row or column of crosstab.

ref：a text corresponds to label or code of arbitrary reference group or a number corresponds to row of group in output (when we want change reference group of output). Also user can enter an arbitrary numeric vector as reference group. Default is Null that used for total of all group as reference (special case that equivalent to the Kruskal-Wallis test).

library(Ridit)

ridit(inptc$eft,inptc$group)  #不设置ref，默认参照组是所有组的综合结果显示：三组无统计学差异（chi2=2.780，P=0.249>0.05）。当某组例数较多时可以将其指定为参照值，命令：

ridit(inptc$eft,inptc$group,ref=3)  #参照水平为对照组10%的硫磺软膏，其实对照组例数并不多，将其设为参照水平并不合适，该命令只做演示如果你的数据已经概况呈频数格式，函数需要指定交叉表和分组变量。

library(readxl)

inptc2 <- read_excel("D:/Temp/NPtdata.xlsx",3)  #从excel中导入第3个sheet的数据，sheet3中的数据是已经概括成频数的格式

crosstable<-xtabs(fre~group+eft,data=inptc2)  #从频数表生成交叉表，频数变量写在公式的左侧。要分析的变量应在公式的右侧（即~符号的右方），第一个变量为行变量，第二个变量为列变量，以+作为分隔符。如果是从原始数据表格生成，命令为crosstable<-xtabs(~group+eft,data=inptc)或crosstable<-table(inptc$group,inptc$eft)

ridit(crosstable,1) #对交叉表crosstable进行Ridit分析，其中行为分组变量（如果列为分组变量则g=2）

如果总体有统计学意义，R中的Ridit并未给出直接的两两比较方法（或者是我没找到），或许你可以分别进行两两Ridit比较，然后通过P值校正来进行。P值校正函数p.adjust {stats}：p.adjust(p, method = p.adjust.methods, n = length(p)) # p.adjust.methods：c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none")

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：多个独立样本的非参数检验及其多重比较。

2020.06.21，夏至

END