《R语言实战》自学笔记69-重抽样和自助法

Python07

《R语言实战》自学笔记69-重抽样和自助法,第1张

数据准备

许多实际情况中统计假设(假定观测数据抽样自正态分布或者其他性质较好的理论分布)并不一定满足,比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况,这时基于随机化和重抽样的统计方法就可派上用场。

置换检验的定义

置换检验(Permutation test),也称随机化检验或重随机化检验,是Fisher于20世纪30年代提出的一种基于大量计算(computationally intensive),利用样本数据的全(或随机)排列,进行统计推断的方法,因其对总体分布自由,应用较为广泛,特别适用于总体分布未知的小样本资料,以及某些难以用常规方法分析资料的假设检验问题。

置换检验的原理

1、提出原假设,比如XX处理后结果没有变化

2、计算统计量,如两组的均值之差,记作t0

3、将所有样本放在一起,然后随机排序进行分组,再计算其统计量t1

4、重复第3步骤,直至所有排序可能性都齐全(比如有A组有n样本,B组有m样本,则总重复次数相当于从n+m中随机抽取n个的次数),得到一系列的统计量(t1-tn)

5、最后将这些统计量按照从小到大排序,构成抽样分布,再看t0是否落在分布的置信区间内(如95%置信区间),这时候可计算一个P值(如果抽样总体1000次统计量中大于t0的有10个,则估计的P值为10/1000=0.01),落在置信区间外则拒绝原假设

6、如果第3步骤是将所有可能性都计算了的话,则是精确检验;如果只取了计算了部分组合,则是近似结果,这时一般用蒙特卡罗模拟(Monte Carlo simulation)的方法进行置换检验

7、置换检验和参数检验都计算了统计量,但是前者是跟置换观测数据后获得的经验分布进行比较,后者则是跟理论分布进行比较。

请牢记:置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样(当做近似检验时)。因此,每次检验的结果都有所不同。

coin包提供了一个进行置换检验的一般性框架。通过该包,你可以回答如下问题。

响应值与组的分配独立吗?

两个数值变量独立吗?

两个类别型变量独立吗?

表12-2列出来的每个函数都是如下形式:

function_name(formula, data, distribution=)

其中:

 formula描述的是要检验变量间的关系。示例可参见表12-2;

 data是一个数据框;

 distribution指定经验分布在零假设条件下的形式,可能值有exact,asymptotic和

approximate。

若distribution = "exact",那么在零假设条件下,分布的计算是精确的(即依据所有可能的排列组合)。当然,也可以根据它的渐进分布(distribution = "asymptotic")或蒙特卡洛重抽样(distribution = "approxiamate(B = #)")来做近似计算,其中#指所需重复的次数。

distribution = "exact"当前仅可用于两样本问题。

传统t检验表明存在显著性差异(p <0.05),而精确检验却表明差异并不显著(p >0.072)。

第7章我用自己的数据进行了t检验,对比一下传统t检验和置换检验,结果如下:

两种检验方式下结果都是显著的

Wilcoxon-Mann-Whitney U检验

coin包规定所有的类别型变量都必须以因子形式编码。

wilcox.test()默认计算的也是精确分布。

K样本检验的置换检验

通过chisq_test()或cmh_test()函数,我们可用置换检验判断两类别型变量的独立性。 当数据可根据第三个类别型变量进行分层时,需要使用后一个函数。若变量都是有序型,可使用lbl_test()函数来检验是否存在线性趋势。

卡方独立性检验

卡方独立性检验的置换检验

你可能会有疑问,为什么需要把变量Improved从一个有序因子变成一个分类因子?(好问题!)这是因为,如果你用有序因子,coin()将会生成一个线性与线性趋势检验,而不是卡方检验。

结果解读:两种检验下p值都是小于0.05,说明Treatment和Improved之间相互不独立

自己数据的演示

结果解读:p值均为1,表明nitrogen和variety相互独立。

spearman_test()函数提供了两数值变量的独立性置换检验。

当处于不同组的观测已经被分配得当,或者使用了重复测量时,样本相关检验便可派上用场。

对于两配对组的置换检验,可使用wilcoxsign_test()函数;多于两组时,使用friedman_test()函数。

自己数据演示

lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版,能够进行置换检验,而非正态理论检验。

lmp()和aovp()函数的参数与lm()和aov()函数类似,只额外添加了perm =参数。

perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样,直至估计的标准差在估计的p值0.1之下,判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意,若观测数大于10,perm = "Exact"将自动默认转为perm = "Prob",因为精确检验只适用于小样本问题。

简单线性回归的置换检验

R语言实战的例子:

多项式回归的置换检验

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

当两种方法所得结果不一致时,你需要更加谨慎地审视数据,这很可能是因为违反了正态性假设或者存在离群点。

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

R语言实战的例子:

自己数据集的例子:

值得注意的是,当将aovp()应用到方差分析设计中时,它默认使用唯一平方和法(SAS也称为类型III平方和)。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和(SAS是类型I平方和)。每种效应依据模型中先出现的效应做相应调整。对于平衡设计,两种方法结果相同,但是对于每个单元格观测数不同的不平衡设计,两种方法结果则不同。不平衡性越大,结果分歧越大。若在aovp()函数中设定seqs = TRUE,可以生成你想要的序贯平方和。

你可能已经注意到,基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好,两种方法结果的一致性也验证了正态理论方法适用于上述示例。

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

置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。不过,置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是,这正是自助法大显神通的地方。

所谓自助法,即从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。 无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。

倘若你假设均值的样本分布不是正态分布,该怎么办呢?可使用自助法。

(1)从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一直都不会被选中。

(2)计算并记录样本均值。

(3)重复1和2一千次。

(4)将1000个样本均值从小到大排序。

(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限定了95%的置信区间。

boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。

一般来说,自助法有三个主要步骤。

(1)写一个能返回待研究统计量值的函数。如果只有单个统计量(如中位数),函数应该返回一个数值;如果有一列统计量(如一列回归系数),函数应该返回一个向量。

(2)为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理。

(3)使用boot.ci()函数获取第(2)步生成的统计量的置信区间。

主要的自助法函数是boot(),它的格式为:

bootobject <- boot(data=, statistic=, R=, ...)

参数见下表:

boot()函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在bootobject中。

你可以用bootobject t0和bootobject t来获取这些元素。

一旦生成了自助样本,可通过print()和plot()来检查结果。如果结果看起来还算合理, 使用boot.ci()函数获取统计量的置信区间。格式如下:

boot.ci(bootobject, conf=, type= )

type参数设定了获取置信区间的方法。perc方法(分位数)展示的是样本均值,bca将根据偏差对区间做简单调整。

回归的R平方值

1000次自助抽样

输出结果

结果可视化

95%的置信区间获取

回归系数向量函数

自助抽样1000次

获得车重和发动机排量95%的置信区间

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

参考资料:

数据准备

函数:legend(location, title, legend, ....)

参数详解:

x和y:用于定位图例,也可用关键词"bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" 和 "center";当图例用关键词设置位置后,inset = 分数,可以设置其相对位置;

legend:指定图例标签,字符或表达式向量;

fill:用特定的颜色进行填充;

col:设置图例中出现的点或线的颜色;

border:当fill = 参数存在的情况下,用于指定填充的边框颜色;

lty, lwd:图例中线的类型与宽度;

pch:点的类型;

angle:阴影的角度;

density:阴影线的密度;

cex:指定图例显示大小;

bg:指定图例的背景色;

bty:指定图例框是否画出,默认o为画出,n为不画出;

box.lty, box.lwd, box.col: 设置图例边框线型,线粗,颜色,box.lty为虚线,box.lwd决定粗线,box.col决定颜色;

pt.bg:图例中点的背景色;

pt.cex:图例中点的大小;

pt.lwd:图例中点边缘的线宽;

x.intersp:图例中文字离图片的水平距离;

y.intersp:图例中文字离图片的垂直距离;

adj:图例中字体的相对位置;

text.width:图例中字体所占的宽度,调整后图例整个宽度也跟着变化了;

text.col:图例字体的颜色;

text.font:图例字体;

merge:逻辑值,merge=TRUE,合并点与线,但不填充图例框,默认为TRUE;

trace:逻辑值,trace=TRUE显示图例信息;

plot:逻辑值,plot=FALSE不画出图例;

ncol:图例中分类的列数;

horiz:逻辑值,horiz=TRUE,水平放置图例;

title:给图例加标题;

xpd:xpd=FALSE,即不允许在作图区域外作图,改为TRUE即可,与par()参数配合使用;

title.col:标题颜色;

title.adj:图例标题的相对位置,0.5为默认,在中间。0最左,1为最右;

seg.len:指定图例中线的线长,长度单位为字符宽度。

1 图例方位

2 修饰图例

3 图例绘制在图外

4 自定义图例

有时候绘制出的图是分组图,这时候需要自定义绘制图例。

参考资料:

《R语言实战》(中文版),人民邮电出版社,2013.

R语言绘制图例(legend)的各种问题_详细综合解析, https://blog.csdn.net/xiangyong58/article/details/54579293

投必得R语言教程,第二讲 R作图-基础-图形参数设置:标题、图例、文字, https://mp.weixin.qq.com/s?__biz=MzU1Mzc3OTIwNg==&mid=2247495531&idx=1&sn=cdd80d4e950ae2b344cf188c68922fa8&chksm=fbef0602cc988f14e93b71bc6fcc2fda782d3ae9a3a67601a6c87756f4ae85bcc5d9c56d9b51&scene=21#wechat_redirect

一幅图解决R语言绘制图例的各种问题, https://blog.csdn.net/weixin_30469895/article/details/96649305

数据准备

t检验,亦称student t检验(Student's t test),主要用于样本含量较小(例如n <30),总体标准差σ未知的正态分布。t检验是用t分布理论来推论差异发生的概率,从而比较两个平均数的差异是否显著。

t检验的适用条件为样本分布符合正态分布。

R中检验正态分布的函数:

shapiro.test()

结果p值要是小于0.05,样本分布是非正态分布,如果大于0.05,样本分布是正态分布。

t检验可分为单总体检验和双总体检验,以及配对样本检验。

单总体t检验是检验一个样本平均数与一个已知的总体平均数的差异是否显著。

个人理解的应用实例:已知一个玉米品种的产量是8000 kg/ha,在一个田间试验中测定这个玉米品种的产量,单样本t检验要做的就是检验田间试验测定的产量与已知产量是否相等。

单样本t检验的假设:

H0:样本均值与已知的总体均值相等。

H1:样本均值与已知的总体均值不相等。

t统计量的计算:

m:样本平均值;

:已知总体的均值;

S:样本标准差,自由度df=n-1。

n:样本量。

单样本t检验R调用函数:

t.test(x, mu, alternative = "two.sided")

x:数据向量;

mu:理论平均值。默认为0,可根据自己统计计算需求更改;

alternative:备择假设。允许值为“two.sided”(默认),也可以根据需要设置为“greater”或“less”之一。

结果解释:p值小于0.05,结论是v1的平均值与理论值1.5有显著差异。

检验两个样本平均数与其各自所代表的总体的差异是否显著。

个人理解的应用实例:检验两个玉米品种产量是否存在差异。

t.test(y ~ x, data)

其中的y是一个数值型变量,x是一个二分变量。

t.test(y1, y2)

其中的y1和y2为数值型向量(即各组的结果变量)。可选参数data的取值为一个包含了这些变量的矩阵或数据框。

t检验默认假定方差不相等,并使用Welsh的修正自由度。你可以添加一个参数var.equal=TRUE以假定方差相等,并使用合并方差估计。默认的备择假设是双侧的(即均值不相等,但大小的方向不确定)。你可以添加一个参数alternative="less"或alternative="greater"来进行有方向的检验。

结果解读:得到结果中P值小于0.05,说明要拒绝原假设(两品种v1值无差异),接受备择假设,即两品种v1值差异显著。

非独立样本的t检验假定组间的差异呈正态分布。

个人理解的应用实例:一个玉米品种接受两个施氮处理,两个施氮处理下玉米的产量是否存在差异。

t.test(y1, y2, paired=TRUE)

其中的y1和y2为两个非独立组的数值向量。

结果解读:不同氮素水平的比较显示p值小于0.05,说明v1值在两个氮水平间差异显著;而两个年份下v1值无显著差异。

如果想在多于两个的组之间进行比较,应该怎么做?如果能够假设数据是从正态总体中独立抽样而得的,那么你可以使用方差分析(ANOVA)。ANOVA是一套覆盖了许多实验设计和准实验设计的综合方法。

参考资料: