【原】R语言统计系列第10篇-自助法

医科研 2021-01-25

展开全文

对于正态分布或其他已知分布的数据，有相应的假设检验与置信区间的计算方法，但是当数据抽样自未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况，就需要使用基于随机化和重抽样的统计方法。

上次推文主要我们介绍了置换检验，本次推文主要介绍自助法。

自助法

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

简单来讲，自助法就是不依赖变量具体分布形式计算置信区间的方法。

自助法步骤：（以下引自R语言实战）

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

（2）计算并记录样本均值；

（3）重复1和2一千次；

（4）将1000个样本均值从小到大排序；

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

针对不服从正态分布的样本均值，自助法的优势十分明显，因此常用于潜在分布未知、出现离群值、样本量过小、没有可选的参数方法等情况。

对单个统计量使用自助法

rsq<-function(formula,data,indices){
  d<-data[indices,]
  fit<-lm(formula,data=d)
  return(summary(fit)$r.square)
}
library(boot)
set.seed(1234)
results<-boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)
print(results)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = mtcars, statistic = rsq, R = 1000, formula = mpg ~ 
##     wt + disp)
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 0.7809306 0.01379126  0.05113904
plot(results)

img

#计算95%的置信区间
boot.ci(results,type=c("perc","bca"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, type = c("perc", "bca"))
## 
## Intervals : 
## Level     Percentile            BCa          
## 95%   ( 0.6753,  0.8835 )   ( 0.6344,  0.8561 )  
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable

我们来分析一下这段代码：

首先，我们定义了一个rsq函数，这个函数表示对于formula计算R方值。此处主要需要修改的代码是 return(summary(fit)$r.square)

其次，我们进行自主法，重复1000次。此处主要需要修改的代码是data和formula

最后，计算95%置信区间，使用“perc”和“bca”的方法。

从图中，我们可以看出自助的R方不呈正态分布。

依据分位数法计算置信区间：Percentile( 0.6753, 0.8835 )

依据偏差对置信区间进行调整：BCa ( 0.6344, 0.8561 )

对多个统计量使用自助法

bs<-function(formula,data,indices){
  d<-data[indices,]
  fit<-lm(formula,data=d)
  return(coef(fit))  
}
library(boot)
set.seed(1234)
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  4.715497e-02 2.546106756
## t2* -3.35082533 -4.908125e-02 1.154800744
## t3* -0.01772474  6.230927e-05 0.008518022
boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)
## 
## 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.2734867271 2.523874212
## t2* -3.35082533 -0.1145055892 1.171195882
## t3* -0.01772474  0.0002178158 0.008599433
plot(results,index=2)

img

## 计算车重回归系数的置信区间
boot.ci(results,type="bca",index=2)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, type = "bca", index = 2)
## 
## Intervals : 
## Level       BCa          
## 95%   (-5.477, -0.937 )  
## Calculations and Intervals on Original Scale
## 计算发动机排量回归系数的置信区间
boot.ci(results,type="bca",index=3)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, type = "bca", index = 3)
## 
## Intervals : 
## Level       BCa          
## 95%   (-0.0334, -0.0011 )  
## Calculations and Intervals on Original Scale