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)的差异有无统计学意义?
多个独立样本,连续型资料,可以考虑的是方差分析。方差分析需要满足正态分布和方差齐性。
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)
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来进行演示。
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)
在Nemenyi法中,Only for method = "Chisq" a tie correction is employed.
posthoc.kruskal.nemenyi.test(inpt$cortisol,inpt$G,dist="Chisquare")
②Dunn法:posthoc.kruskal.dunn.test {PMCMR}
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, ...)
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"
示例2
244例疥疮患者接受2.5%和1%的二硫化硒洗剂、10%的硫磺软膏治疗,其中有效药物硫磺软膏外用作为对照,考察二硫化硒治疗疥疮的疗效。
岭南皮肤性病科杂志,2005,12(1):57-58.
在很多杂志期刊中,此类数据常误用卡方来进行分析。这种分组变量为有序/无序分类变量、指标变量为有序分类变量的资料,应该使用非参数检验法。三组独立样本,宜采用前面刚刚讲过的Kruskal-Wallis检验。文中使用的并不是秩和检验法,而是Ridit法,我们也会在此演示一下。
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)
library(PMCMR)
posthoc.kruskal.dunn.test(eft~group,data=inptc)
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,默认参照组是所有组的综合
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")
原创不易,欢迎“在看”!
2020.06.21