R语言置换检验

R 专栏收录该内容
3 篇文章 0 订阅

原文地址:http://blog.csdn.net/gdyflxw/article/details/54141629

置换检验

  双样本均值检验的时候,假设检验的方法就是,检查正态性、独立性、方差齐性,分别对应的参数非参数方法进行假设检验,但是,这些方法都要求样本数必须有多少多少,但是,由于试验时,各种条件的限制,导致样本量过小,此时以上方法几乎都会失真,置换检验就应运而生了。 
  Permutation test 置换检验是Fisher于20世纪30年代提出的一种基于大量计算 (computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。在具体使用上它和Bootstrap Methods类似,通过对样本进行顺序上的置换,重新计算统计检验量,构造经验分布,然后在此基础上求出P-value进行推断。 
  置换检验的操作方法:假设有两组待检数据,A组有m个数据,B组有n个数据,均值差为d0,现把所有数据放在一起进行随机抽取,抽出m个放入A组,剩下n个放入B组,计算A、B两组的均值差记为d1,再放在一起进行随机重抽m、n两组,得到均值差记为d2,重复这个步骤k次得到(d3……dk),于是d1……dk可以画出一张正态图,然后看d0落在什么方,若落在置信水平之外,即可以显著说明它们是有差异的。 
  R代码如下:

a<-c(24,43,58,67,61,44,67,49,59,52,62,50,42,43,65,26,33,41,19,54,42,20,17,60,37,42,55,28)
group<-factor(c(rep("A",12),rep("B",16)))
data<-data.frame(group,a)
find.mean<-function(x){
    mean(x[group=="A",2])-mean(x[group=="B",2]) 
} 
results<-replicate(999,find.mean(data.frame(group,sample(data[,2])))) 
p.value<-length(results[results>mean(data[group=="A",2])-mean(data[group=="B",2])])/1000
hist(results,breaks=20,prob=TRUE)
lines(density(results))

   
From:https://www.plob.org/article/3176.html

coin包置换检验

coin包介绍

  coin包中的置换检验有以下几种:

检 验 coin函数
两样本和K样本置换检验 oneway_test(y ~ A)
含一个分层(区组)因子的两样本和K样本置换检验 oneway_test(y ~ A | C)
Wilcoxon-Mann-Whitney秩和检验 wilcox_test(y ~ A)
Kruskal-Wallis检验 kruskal_test(y ~ A)
Person卡方检验 chisq_test(A ~ B)
Cochran-Mantel-Haenszel检验 cmh_test(A ~ B | C)
线性关联检验 lbl_test(D ~ E)
Spearman检验 spearman_test(y ~ x)
Friedman检验 friedman_test(y ~ A | C)
Wilcoxon符号秩检验 wilcoxsign_test(y1 ~ y2)

注:在上表中,y和x是数值变量,A和B是分类因子,C是类别型区组变量,D和E是有序因子,y1和y2是相匹配的值变量 
表中所有的函数使用方法都一样:

functionName(formula,dataframe,distribution),其中distribution指定经验分布在零假设条件下的形式,可能值有exact,asymptotic和approximate,若distribution = "exact",那么在零假设条件下,分布的计算是精确的(即依据所有可能的排列组合)。当然,也可以根据它的渐进分布(distribution = "asymptotic")或蒙特卡洛重抽样(distribution = "approxiamate(B = #)")来做近似计算,其中#指所需重复的次数。distribution = "exact"当前仅可用于两样本问题。

原函数与置换检验比较

函数 简介 程序及结果
t.test() 双样本均值t检验 > score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)
> treatment <- factor(c(rep(“A”, 5), rep(“B”, 5)))
> mydata <- data.frame(treatment, score)
> t.test(score ~ treatment, data = mydata, var.equal = TRUE)
          Two Sample t-test
data: score by treatment
t = -2.345, df = 8, p-value = 0.04705
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  -19.0405455    -0.1594545
sample estimates:
mean in group A mean in group B
     51.0     60.6
oneway_test() 双样本均值置换检验 > oneway_test(score ~ treatment, data = mydata, distribution = “exact”)
    Exact Two-Sample Fisher-Pitman Permutation Test
data: score by treatment (A, B)
Z = -1.9147, p-value = 0.07143
alternative hypothesis: true mu is not equal to 0
wilcox.test() 双样本秩和独立性检验 > wilcox.test(Prob~So,data=UScrime)
     Wilcoxon rank sum test
data: Prob by So
W = 81, p-value = 8.488e-05
alternative hypothesis: true location shift is not equal to 0
wilcox_test() 双样本秩和独立性置换检验 > UScrime2 <- transform(UScrime, So = factor(So))
> wilcox_test(Prob ~ So, data = UScrime2, distribution = “exact”)
    Exact Wilcoxon-Mann-Whitney Test
data: Prob by So (0, 1)
Z = -3.7493, p-value = 8.488e-05
alternative hypothesis: true mu is not equal to 0
aov() 单因素方差分析 > library(multcomp)
>summary(aov(response~trt,data=cholesterol))
  Df Sum Sq  Mean Sq  F value Pr(>F)
trt 4 1351.4   337.8    32.43  9.82e-13 ***
Residuals 45 468.8 10.4
oneway_test() K样本置换检验 > oneway_test(response ~ trt, data = cholesterol, distribution = approximate(B = 9999))
  Approximative K-Sample Fisher-Pitman Permutation Test
data: response by
trt (1time, 2times, 4times, drugD, drugE)
chi-squared = 36.381, p-value < 2.2e-16
chisq.test() 卡方列联表均值差异检验 > chisq.test(xtabs(~Treatment+Improved,Arthritis))
   Pearson’s Chi-squared test
data: xtabs(~Treatment + Improved, Arthritis)
X-squared = 13.055, df = 2, p-value = 0.001463
chisq_test() 卡方置换检验 > chisq_test(Treatment ~ Improved, data = transform(Arthritis, Improved = as.factor(as.numeric(Improved))),distribution = approximate(B = 9999))
   Approximative Pearson Chi-Squared Test
data: Treatment by Improved (1, 2, 3)
chi-squared = 13.055, p-value = 0.0012
mantelhaen.test() 分层卡方检验,看是否把相关因素划分出去 > mytable <- xtabs(~Treatment+Improved+Sex, data=vcd::Arthritis)
> mantelhaen.test(mytable)
    Cochran-Mantel-Haenszel test
data: mytable
Cochran-Mantel-Haenszel
M^2 = 14.632, df = 2, p-value = 0.0006647
cmh_test() 分层卡方置换检验,看是否把相关因素划分出去 > cmh_test(mytable)
   Asymptotic Generalized Cochran-Mantel-Haenszel Test
data: Improved by
Treatment (Placebo, Treated) 
stratified by Sex
chi-squared = 14.632, df = 2, p-value = 0.0006647
cor() spearman等级相关系数 > with(states,cor(Illiteracy,Murder,method=”spearman”))
[1] 0.6723592
spearman_test() 数值独立性置换检验(两数值变量独立即不相关) > spearman_test(Murder~Illiteracy,data=states)
   Asymptotic Spearman Correlation Test
data: Murder by Illiteracy
Z = 4.7065, p-value = 2.52e-06
alternative hypothesis: true rho is not equal to 0
t.test(paired=T) 非独立样本的配对t检验,检验均值是否相等 > with(MASS::UScrime,t.test(U1,U2,paired=TRUE))
     Paired t-test
data: U1 and U2
t = 32.407, df = 46, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
57.67003 65.30870
sample estimates:
mean of the differences 
61.48936
wilcoxsign_test() wilcox符号秩置换检验,检验均值是否相等 > wilcoxsign_test(U1 ~ U2, data = MASS::UScrime,distribution = “exact”)
   Exact Wilcoxon-Pratt Signed-Rank Test
data: y by x (pos, neg) 
stratified by block
Z = 5.9691, p-value = 1.421e-14
alternative hypothesis: true mu is not equal to 0
friedman_test() 多组别独立性置换检验,检验均值是否相等 > USc<-MASS::UScrime[,c(“U1”,”U2”)]
> USc$U3<-sample(as.matrix(USc),47)
>friedman_test(value~variable|ID,data=transform(reshape::melt(data.frame(USc,ID=seq(1,47)),id.vars=”ID”),ID=as.factor(ID)))
      Asymptotic Friedman Test
data: value by
variable (U1, U2, U3) 
stratified by ID
chi-squared = 51.384, df = 2, p-value = 6.953e-12

  coin包的介绍至此结束,当然还有一个lbl_test()函数未列出,暂时还不晓得有什么用,以后再说。

lmPerm包置换检验

lmPerm包介绍

  lmPerm包可以做非正态理论检验,包含的函数为lmp()以及aovp()两个,它们与lm()和aov()类似,只是多了一个perm参数(perm=”Exact”,”Prob”,”SPR”),参数值”Exact”根据所有可能的排列组合生成精确检验,”Prob”从所有可能的排列中不断抽样,直至估计的标准差在估计的p值0.1之下,判停准则由可选的Ca参数控制,SPR使用贯序概率比检验来判断何时停止抽样。若观测数大于10,perm=”Exact”会自动转化为perm=”Prob”,因为精确检验只适用于小样本问题。 
  因为只涉及了两个函数,这个包就不贴代码和结果,仅说明一下差异是什么,

回归(简单、多项式、多元)

  首先是lm与lmp,除了函数的用法多了个perm参数之外,所得结果模板(注意,是模板,而非结果,结果出现差异应该去找数据的问题,如两者结果不一致,则需要重新审视数据的可靠性)存在差异: 
  1)少了常数项,但可以通过各变量均值求得,注意,使用coefficients(fit)所得的常数项是错的! 根据回归线必过均值点的定义,可以使用各变量的均值来计算其常数项。如多元分析中的例子计算方式为:

mean(states$Murder)-sum(colMeans(states)[names(coefficients(fit)[c(-1)])]*(coefficients(fit)[c(-1)]))

  2)回归系数项中多了Iter一栏,它表示要达到判停准则所需要的迭代次数。

方差分析

  与回归一致,所有使用aov分析的地方都可以使用aovp来代替,区别就是,aov用的是F统计量,而aovp使用的是置换法,Iter为判停准则的迭代次数。 
  需要注意的是,aovp使用的是唯一平方和方法,每种效应根据其它效应进行调整,而aov使用的是序贯平方平法,每种效应根据先出现的效应进行调整,这两个方法在不平衡设计中所得结果不同,越不平衡的设计,差异越大。可以在aovp函数里加入参数seqs=TRUE可以生成序贯平方和的计算结果。 
  

点评

  置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。 
  

自助法

  置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。不过,置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是,这正是自助法大显神通的地方。 
  自助法的步骤: 
  1. 一个样本数为n的样本,进行m次有放回抽样; 
  2. 计算并记录样本统计量(比如均值、方差、甚至t检验量等,可以一个,可以多个); 
  3. 重复1000到2000次,或者更多,并把它们从小到大进行排序; 
  4. 根据双尾95%分位点,即2.5%和97.5%分位数,即为95%置信区间的下限和上限。

boot包

  boot包可以进行自助法抽检,并生成相应的置信区间。 
  主要的步骤如下: 
  1. 定义函数,返回一个统计值或一个向量(多个统计值),函数要包括indices参数,以便boot()函数用它从每个重复中选择实例,主要是stype参数,默认为i(索引值),还有f(频率)和w(权重),indices可以简定为i; 
  2. 用boot(data,sitisctic,R,……)函数生成一个bootobject。 
  3. 使用boot.ci(bootobject,conf,type)生成置信区间,其中conf定义置信区间,type定义置信区间类型(即计算方法),包含norm、basic、stud、perc、bca和all(其中norm为正态分布的置信区间计算方法,约两个标准差距离,perc为上下分位数计算方法,stud为t分布计算方法),若返回值为向量,则利用index参数来指定某个变量的置信区间。 
  4. 其它相关数据:比如bootobject t0bootobject t为重复R次的统计量值(一个“R*统计量个数”的矩阵)

  最后谨记:置换检验和自助法并不是万能的,它们无法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳,或者样本量过小而无法准确地反映总体情况,这些方法也是爱莫能助。

  • 2
    点赞
  • 1
    评论
  • 14
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值