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

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

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

自助法

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

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

自助法步骤：（以下引自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


本例计算的是截距、车重和发动机排量的回归系数的置信区间。

注意事项：

①初始样本大小10~30即可得到足够好的结果。

②重复1000次一般可满足需求。

注意，如果样本代表性不佳，或者样本量过小无法反映总体情况，使用基于随机化和重抽样的统计方法也无法将其转化为好数据，因此分析数据的准确是结论准确的前提。

好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！