12.4 置换检验点评
除coin和lmPerm包外,R还提供了其他可做置换检验的包。perm包能实现coin包中的部分功能,因此可作为coin包所得结果的验证。corrperm包提供了有重复测量的相关性的置换检验。
logregperm包提供了Logistic回归的置换检验。另外一个非常重要的包是glmperm,它涵盖了广义线性模型的置换检验依靠基础的抽样分布理论知识,置换检验提供了另外一个十分强大的可选检验思路。对于上面描述的每一种置换检验,我们完全可以在做统计假设检验时不理会正态分布、t分布、F分布或者卡方分布。当然,置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。
12.5 自助法
所谓自助法,即从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。
无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。比如,你想计算一个样本均值95%的置信区间。假设均值的样本分布不是正态分布:
(1) 从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一
直都不会被选中。
(2) 计算并记录样本均值。
(3) 重复1和2一千次。
(4) 将1000个样本均值从小到大排序。
(5) 找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限
定了95%的置信区间。
12.6 boot 包中的自助法
boot包扩展了自助法和重抽样的相关用途。可以对一个统计量(如中位数)或一个统计
量向量(如一列回归系数)使用自助法.
自助法有三个主要步骤。
(1) 写一个能返回待研究统计量值的函数。如果只有单个统计量(如中位数),函数应该返回
一个数值;如果有一列统计量(如一列回归系数),函数应该返回一个向量。
(2) 为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理。
(3) 使用boot.ci()函数获取第(2)步生成的统计量的置信区间。
主要的自助法函数是boot(),它的格式为:bootobject<-boot(data=,statistic=,R=,…)
data:量、矩阵或者数据框
statistic:生成k个统计量以供自举的函数(k=1时对单个统计量进行自助抽样)函数需包括indices参数,以便boot()函数用它从每个重复中选择实例
R:自助抽样的次数
...:其他对生成待研究统计量有用的参数,可在函数中传输
boot()函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指
标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在
bootobject中。boot()函数中返回对象所含的元素
t0 从原始数据得到的k个统计量的观测值
t 一个R × k矩阵,每行即k个统计量的自助重复值
你可以如bootobject$t0和bootobject$t这样来获取这些元素。
一旦生成了自助样本,可通过print()和plot()来检查结果。如果结果看起来还算合理,
使用boot.ci()函数获取统计量的置信区间。格式如下:
Boot.ci(bootobject,couf=type=)
bootobject boot()函数返回的对象
conf 预期的置信区间(默认:conf =0.95)
type 返回的置信区间类型。可能值为norm、basic、stud、perc、bca和all(默认:type =all)
type参数设定了获取置信区间的方法。perc方法(分位数)展示的是样本均值,bca将根据
偏差对区间做简单调整
12.6.1 对单个统计量使用自助法
> rsq<-function(formula,data,indices){
+ d<-data[indices,]
+ fit<-lm(formula,data=d)
+ return(summary(fit)$r.square)
+ }
> set.seed(1234)
> results<-boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)
> print(results)
ORDINARY NONPARAMETRICBOOTSTRAP
Call:
boot(data = mtcars,statistic = rsq, R = 1000, formula = mpg ~
wt + disp)
Bootstrap Statistics :
original bias std. error
t1* 0.78093060.0133367 0.05068926
可以看到,自助的R平方值不呈正态分布。它的95%的置信区间可以通过如下代
码获得:
> boot.ci(results,type=c("perc","bca"))
BOOTSTRAP CONFIDENCE INTERVALCALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot.ci(boot.out = results, type =c("perc", "bca"))
Intervals :
Level Percentile BCa
95% ( 0.6838, 0.8833 ) ( 0.6344, 0.8549 )
Calculations and Intervals on OriginalScale
Some BCa intervals may be unstable
12.6.2 多个统计量的自助法
首先,创建一个返回回归系数向量的函数:
> bs<-function(formula,data,indices){
+ d<-data[indices,]
+ fit<-lm(formula,data=d)
+ return(coef(fit))
+ }
> results<-boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)
> print(results)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = mtcars, statistic = bs, R= 1000, formula = mpg ~
wt + disp)
Bootstrap Statistics :
original bias std. error
t1* 34.96055404 0.1101475088 2.445950503
t2* -3.35082533 -0.08573079461.114315903
t3* -0.01772474 0.0003223228 0.008201569
plot(results,index=2)
12.7 小结
本章,我们介绍了一系列基于随机化和重抽样的计算机密集型方法,它们使你无需理论分布
的知识便能够进行假设检验,获得置信区间。当数据来自未知分布,或者存在严重的离群点,或
者样本量过小,又或者没有参数方法可以回答你感兴趣的假设问题时,这些方法是非常实用的。