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

(0)

相关推荐

  • 别人的膀胱癌ceRNA如何做出新花样?

    Identification of a potentially functional circRNA–miRNA–mRNA regulatory network for investigating p ...

  • Spring Boot 论坛开发项目

    ┃  ┃  ┃  ┣━62 62 项目发布 - 服务器选购.mp46 P! g: N! M! X ┃  ┃  ┃  ┣━01 1 开篇介绍.mp4 ┃  ┃  ┃  ┣━47 47 规范提问的标签库. ...

  • R语言统计系列第13篇-K-M生存曲线与logrank检验

    白介素2的读书笔记,分享临床科研干货,一起见证时间的力量 [科研绘图需求点我][付费精品合集][SEER点我] [临床预测模型专辑点我][生物信息数据库挖专辑点我] [临床数据挖掘专辑点我][SEER ...

  • R语言统计系列第12篇-泊松回归

    [科研绘图点我][付费精品合集][SEER点我] 今天是各类统计方法R语言实现的第12期,我们主要介绍泊松回归. 泊松回归 计数型,如某地区某年发生肿瘤患者的人数等,常用泊松回归,假设因变量Y服从泊松 ...

  • R语言统计系列第11篇-Logistic回归

    [科研绘图点我][付费精品合集][SEER点我] 今天是各类统计方法R语言实现的第11期,我们主要介绍Logistic回归.Logistic回归属于广义线性回归,因此我们从广义线性回归讲起. 广义线性 ...

  • R语言统计系列第9篇-置换检验

    [科研绘图点我][付费精品合集][SEER点我] 对于正态分布或其他已知分布的数据,有相应的假设检验与置信区间的计算方法,但是当数据抽样自未知或混合分布.样本量过小.存在离群点.基于理论分布设计合适的 ...

  • 【R语言统计合集】R语言统计系列10篇推文汇总目录

    [科研绘图点我][付费精品合集][SEER点我] 转眼间,我们的R语言统计系列已经经历了10期,我们从统计学最基础的统计描述讲到置换检验和自主法等内容,算是完成了统计的入门. 在这篇推文里,我一方面将 ...

  • 望闻问切系列第10篇:科学语言解读洪脉、微脉、紧脉、缓脉的机理

    后三期发文预告:病脉所揭示的疾病现象和身体的状态(4)-(6). 发文规则:周一.周三.周五更新.点击关注@三甲医院中医袁大夫查看更多文章. 望闻问切系列文章收录于书籍<深入浅出方解伤寒> ...

  • R语言学习系列之“多变的热图”

    咱公众号也不能只做一个系列,所以经过深思熟虑,打算将来慢慢增加一些内容,主要有以下几个系列TCGA数据分析系列GEO数据分析系列"老板给一个基因,我该怎么办"系列文献阅读系列R语言 ...

  • R语言学习系列之separate {tidyr}

    separate {tidyr}的功能是使用正则表达式分隔符将字符列分隔成多个列.给定正则表达式或字符位置向量,separate()将单个字符列转换为多个列.这个功能在生信分析中用到的频率太多了,比如 ...

  • 望闻问切系列第7篇:现代科学语言解读脉法的基本原理

    后三期发文预告:病脉所揭示的疾病现象和身体的状态(1).(2).(3). 发文规则:周一.周三.周五更新.点击关注@三甲医院中医袁大夫查看更多文章. 望闻问切系列文章收录于书籍<深入浅出方解伤寒 ...