R语言入门--第十一节(置换检验与自助法求置信区间)

Python040

R语言入门--第十一节(置换检验与自助法求置信区间),第1张

原理参考 文章 ,主要思想我认为是求出所有分布的可能(中间的一般为零假设),出现这种分布的概率。

distribution= 参数可为exact(精确模式,即依据所有可能的排列组合,仅适用于两样本问题)、approxiamate(nresample=#)(蒙特卡洛抽样,#指需要重复的次数)、asymptotic(渐进分布抽样)

lmPerm包更擅长方差分析。示例实验设计仍为5组接受不同治疗方法的多组结果比较。

实验示例仍为关节炎的治疗(两种)与效果(无、部分、显著)间的关系

实验示例为研究文盲率与谋杀率是否相关

主要为 lmp() 、 aovp() 两个函数分别对应参数法的 lm() 线性回归、 aov() 方差分析。主要格式上的区别是添加了 perm= 参数。可以为Exact(精确模式)、Prob(不断从可能的序列中抽样,直至估计的标准差在估计的p值0.1之下)、SPR(使用贯序概率比检验来判断何时停止抽样)。值得注意的是当样本观测大于10,perm="Exact"自动默认转为"Prob",因此精确检验只适用于小样本问题。

(1)简单线性回归

实验示例仍为以身高预测体重的设计

(2)多项式回归

高精度拟合身高体重回归关系

(3)多元回归

探究谋杀率与多因素的回归关系

(1)单因素方差分析

(2)单因素协方差分析

实验示例仍为药物对刚出生小鼠体重影响,协变量为怀孕时间

(3)双因素方差分析(交互效应)

实验示例:两种药物分别在不同剂量下对小鼠牙齿长度的影响。

核心思想是有放回的抽样多次(1000次)

(1)写一个能返回带研究统计量的函数;

(2)确定重复数,使用 boot() 函数处理;(一般重复1000次即可;此外有人认为初始样本大小为20-30即可得到足够好的结果);

(3) boot.ci() 函数计算统计量置信区间。

实验示例:使用mtcar数据框,采用多元回归,根据车重和发动机排量来预测汽车的每加仑行驶的英里数。想获得95%的R平方值(预测变量对响应变量可解释的方差比)的置信区间

(1)首先写函数

(2)然后使用boot()函数

(3)最后boot.ci()函数求置信区间

实验示例:使用mtcar数据框,采用多元回归,根据车总和发动机排量来预测汽车的每加仑行驶的英里数。想获取一个统计量向量--三个回归系数(截距项、车总、发动机排量)95%的置信区间。

Permutation test 可以称作是置换检验,Fisher于20世纪30年代提出的一种基于大量计算(computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。在具体使用上它和Bootstrap Methods类似,通过对样本进行顺序上的置换,重新计算统计检验量,构造经验分布,然后在此基础上求出P-value进行推断。

总体上其实就是重新采样增加小样本的整体样本量,然后看其的概率分布来预测假设是否成立。这个基于t-test的,类似于与t-test,但是我觉的在样本量非常少的情况下,用置换检验可以更好的说明问题。这个p-value主要是sobs值在置换经验中均值的分布概率,最终也是利用p-value的值来判断假设是否成立的,看p值的大小,p值小于0.05时,是说明拒绝H0,大于0.05.则是说明服从0假设。

下面说一下,置换检验在R中应用。

排列测试在实验研究中特别相关,我们常常对治疗组之间无差异的拒绝零假设感兴趣。在这些情况下,置换检验很好的地代表了我们的推理过程,因为我们的零假设是两个治疗组在结果上没有差异(即,结果是独立于治疗分配而观察到的)。当我们在测试期间置换结果值时,因此我们看到我们可能具有的所有可能的替代治疗分配排列以及我们观察到的数据的平均差异(相对于我们可以看到的结果是独立的所有差异的治疗任务的位置)。虽然排列测试要求我们看到数据的所有可能排列(可能变得非常大),通过简单地进行大量的重采样,我们可以轻松地进行“近似置换测试”。在期望中,该过程应该近似于排列分布。假设我们的研究中有20个单位那么他的排列的数量是:

这个时候,数字完全超出了我们合理计算的数字,但是我们可以从该排列分布中随机抽样以获得近似排列分布,只需运行大量重新采样即可。让我们看一下使用一些组成数据的例子:

当我们重新取样而无需替换再次重新计算差异:为-0.2612

在这里,我们使用置换处理向量s计算差异并找到非常小的差异而不是用tr。如果我们重复这个过程很多次,我们可以建立我们的近似置换分布(即均值差的采样分布)。我们将使用replicatedo重复我们的排列过程。结果将是每个排列(即我们的分布)的差异向量:

然后我们可以用hist函数看一下这个分布然后画一个直方图来看它的差异

然后我们可以量化这个结果,也就是生成一个p-value值,通过p-value来更直观的观测结果:

2:在R中我们不可能每次都要构造自己的置换检验的分布集, R那我们可以这用coin包中的independence_test函数,但是一个问题是实际问题考虑的可能只是单侧的置换检验的结果:下面是这个函数针对于上述例子的应用。

对于这个函数的应用,它的具体的参数及应用如下

independence_test(asat ~ group, data = asat,

                  ## exact null distribution

                  distribution = "exact",

                  ## one-sided test

                  alternative = "greater",

                  ## apply normal scores to asat$asat

                  ytrafo = function(data)

                      trafo(data, numeric_trafo = normal_trafo),

                  ## indicator matrix of 1st level of asat$group

                  xtrafo = function(data)

                      trafo(data, factor_trafo = function(x)

                          matrix(x == levels(x)[1], ncol = 1)))

除此之外,还可以用oneway_test做,效果类似。

但是oneway_test计算的默认是双侧的P值, 这时的计算最好是用Deducer包中的perm.t.test函数

这个函数的可以很简单的设置是单边的还是双边的,简单的例子使用如下

library(Deducer)

下面简单介绍下perm.t.test的 参数

如下所示

perm.t.test(x,y,statistic=c("t","mean"),alternative=c("two.sided", "less", "greater"), midp=TRUE, B=10000)

Arguments

x  第一个参数向量(数字类型)

y   第二个参数向量(数字类型)

statistic 统计用的标准,t或者均值

alternative 进行统计排列的方式,主要有三种,如上面所示

midp 确定p-value是否应用

B    进行随机置换取样的次数

这就是置换检验在R中的简单介绍