《R语言实战》第二版碎碎念之三:重抽样与自助法

重抽样与自助法。数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况,基于随机化和重抽样的统计方法可派上用场。

(1)置换检验,也称为随机化检验或重随机化检验。coin包对于独立性问题提供了一个非常全面的置换检验的框架,而lmPerm包专门用来做方差分析和回归分析的置换检验。lmPerm包需要用RTools安装。

暂时没搞懂。。。

(2)自助法。从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。boot包扩展了自助法和重抽样的相关用途。自助法一般三个步骤:

①写一个能返回待研究统计量值的函数,如果只有单个统计量,函数应该返回一个数值,如果有一列统计量,函数应该返回一个向量。

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

③使用boot.ci()函数获取步骤2生成的统计量的置信区间。格式为Bootobject<-boot(data=,statistic=,R=,…),其中data可以是向量、矩阵或数据框,statistic生成k个统计量以供自举的函数(k=1时对单个统计量进行自助抽样)。函数需包括indices参数,以便boot()函数用它从每个重复中选择实例。R是自助抽样的次数。boot()函数的返回对象bootobject中包含两个元素,t0是从原始数据得到的k个统计量的观测值,t是一个R×k矩阵,每行即k个统计量的自助重复值。使用boot.ci(bootobject,conf=,type=)函数获取统计量的置信区间,conf是预期的置信区间,默认conf=0.95,type返回置信区间类型,可能值为norm,basic,stud,perc,bca和all,默认type="all"。

对单个统计量使用自助法:

以mtcars数据集为例,根据车重和发动机排量预测每加仑行驶的英里数,获得95%的R平方值的置信区间(预测变量对响应变量可解释的方差比)。

写一个获取R平方值的函数:

rsq <- function(formula, data, indices) {

注:indices必须声明,因为boot()要用其来选择样本。

d <- data[indices,]

fit <- lm(formula, data=d)

return(summary(fit)$r.square)

}

函数返回回归的R平方值。 做大量的自助抽样(比如1000),代码如下:

library(boot)

set.seed(1234)

results <- boot(data=mtcars, statistic=rsq,R=1000, formula=mpg~wt+disp)

boot的对象可以输出 print(results),也可用plot(results)来绘制结果。获取95%的置信区间: boot.ci(results, type=c("perc", "bca")),生成置信区间的不同方法将会导致获得不同的区间。如果0在置信区间外,零假设H0: R平方值=0被拒绝。

多个统计量的自助法:

继续该例,让我们获取一个统计量向量——三个回归系数(截距项、车重和发动机排量) ——95%的置信区间。首先,创建一个返回回归系数向量的函数:

bs <- function(formula, data, indices) {

d <- data[indices,]

fit <- lm(formula, data=d)

return(coef(fit))

}

然后使用该函数自助抽样1000次:

library(boot)

set.seed(1234)

results <- boot(data=mtcars, statistic=bs,R=1000, formula=mpg~wt+disp)

查看结果:print(results)

当对多个统计量自助抽样时,添加一个索引参数,指明plot()和boot.ci()函数所分析bootobject$t的列。在本例中,索引1指截距项,索引2指车重,索引3指发动机排量。如下代码即用于绘制车重结果:plot(results, index=2)

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

boot.ci(results, type="bca", index=2)

boot.ci(results, type="bca", index=3)

注意 在先前的例子中,我们每次都对整个样本数据进行重抽样。如果假定预测变量有固定水平(如精心设计的实验),那么我们最好仅对残差项进行重抽样。参考Mooney & Duval(1993, pp.16-17),其中有简单的解释和算法。

发表评论

电子邮件地址不会被公开。 必填项已用*标注