R语⾔实战笔记--第⼗⼆章重抽样（置换检验）与⾃助法
R语⾔实战笔记–第⼗⼆章重抽样（置换检验）与⾃助法
标签（空格分隔）： R语⾔ 重抽样 ⾃助法 置换检验
置换检验
  双样本均值检验的时候，假设检验的⽅法就是，检查正态性、独⽴性、⽅差齐性，分别对应的参数⾮参数⽅法进⾏假设检验，但是，这些⽅法都要求样本数必须有多少多少，但是，由于试验时，各种条件的限制，导致样本量过⼩，此时以上⽅法⼏乎都会失真，置换检验就应运⽽⽣了。
  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)
mean(x[group=="A",2])-mean(x[group=="B",2])
}
results<-replicate(an(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))
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 = "approxiamat e(B = #)"）来做近似计算，其中#指所需重复的次数。distribution = "exact"当前仅可⽤于两样本问题。
原函数与置换检验⽐较
函数简介程序及结果
> 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检验
> 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
和独⽴性
检验
> st(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
表均值差
异检验
> st(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
检验，看
是否把相
关因素划
分出去
> mytable <- xtabs(~Treatment+Improved+Sex, data=vcd::Arthritis)
> st(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)
bootstrap检验方法  　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
检验，检
验均值是
否相等
> with(MASS::st(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. 其它相关数据：⽐如bootobjectt为重复R次的统计量值（⼀个“R*统计量个数”的矩阵）
  最后谨记：置换检验和⾃助法并不是万能的，它们⽆法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳，或者样本量过⼩⽽⽆法准确地反映总体情况，这些⽅法也是爱莫能助。