R语⾔与回归分析学习笔记（bootstrapmethod）
Bootstrap⽅法在之前的博⽂⾥有提到过，简⽽⾔之，bootstrap⽅法就是重抽样。为什么需要bootstrap⽅法呢？因为bootstrap ⽅法使得我们⽆需分布理论的知识也可以进⾏假设检验，获得置信区间。当数据来⾃未知分布，或者存在严重异常点，⼜或者样本量过⼩，没有参数⽅法解决问题时，bootstrap⽅法将是⼀个很棒的⽅法。
对于回归分析⽽⾔，bootstrap⽆疑对回归的正态性假设做了极⼤地放松，使得回归推断越来越好⽤，也更具有说服⼒。
从博⽂⾥可以看到，对于参数统计，特别是在已知分布的参数估计，bootstrap并没有多⼤的意义，它的结果和矩估计或者极⼤似然估计的结果并没有多⼤的差别（如果有差别会令⼈觉得很奇怪，不是吗？）
Boot包中提供了做bootstrap的两个⼗分好⽤的函数：boot（），boot.ci（）。两者的调⽤格式与参数说明如下：
Boot（）函数：
boot(data, statistic, R, sim ="ordinary", stype = c("i", "f", "w"),
strata= rep(1,n), L = NULL, m = 0, weights = NULL,
< = function(d, p) d, mle = NULL, simple = FALSE, ...,
parallel = c("no", "multicore", "snow"),
ncpus = getOption("pus", 1L), cl = NULL)
参数说明：
Data：数据，可以是向量，矩阵，数据框
Statistic：统计量，如均值，中位数，回归参数，回归⾥的R^2等
R：调⽤统计量函数次数
Boot（）的返回值：
T0：从原始数据中得到的k个统计量的观测值
T:⼀个R*K的矩阵
Boot.ci（）函数：
boot.ci(boot.out, conf = 0.95, type = "all",
index = 1:min(2,length(boot.out$t0)), var.t0 = NULL,
var.t = NULL, t0 = NULL, t = NULL, L = NULL,
h = function(t) t, hdot = function(t) rep(1,length(t)),
hinv = function(t) t, ...)
参数说明：
Boot.out（）：boot函数的返回值
Type：返回置信区间的类型，R中提供的有"norm" ,"basic", "stud","perc", "bca"
⼀、    对单个统计量使⽤bootstrap⽅法
我们以R中的数据集women为例说明这个问题。数据集women列出了美国妇⼥的平均⾝⾼和体重。以体重为响应变量，⾝⾼为协变量进⾏回归，获取斜率项的95%置信区间。
R可以通过以下代码告诉我们答案：
library(boot)
beta<-function(formula,data,indices){
d<-data[indices,]
fit<-lm(formula,data=d)
return(fit$coef[2])
}
result<-boot(data=women,statistic=beta,R=500,formula=weight~height)
boot.ci(result)
输出结果：
BOOTSTRAPCONFIDENCE INTERVAL CALCULATIONS
Basedon 500 bootstrap replicates
CALL:
boot.ci(boot.out= result)
Intervals:
Level      Normal              Basic
95%  ( 3.218, 3.686 )  ( 3.231,  3.704 )
Level    Percentile            BCa
95%  ( 3.196, 3.669 )  ( 3.199,  3.675 )
Calculationsand Intervals on Original Scale
他们与传统的估计差别⼤吗？我们来看看传统的区间估计：
confint(lm(weight~height,data=women))
输出结果：
2.5 %    97.5 %
(Intercept) -100.342655 -74.690679
height        3.253112  3.646888
可以看出，差别并不是很⼤，究其原因，⽆外乎正态性得到了很好的满⾜。我们看其qq图：
很清楚也很显然。Shapiro检验也说明了这样⼀个事实：
Shapiro-Wilk normality test
data:  women$weight
W = 0.9604,p-value = 0.6986
我们在来看⼀个差别较⼤的例⼦：
.    考虑R中的数据集faithful。以waiting为响应变量，eruptions为协变量，建⽴简单回归模型y=α+βx+e。考虑β的95%置信区间。      重复上⾯的步骤，R代码如下：
result<-boot(data=faithful,statistic=beta,R=500,formula=waiting~eruptions)
boot.ci(result)
confint(lm(waiting~eruptions,data=faithful))
qqPlot(lm(waiting~eruptions,data=faithful))
输出结果：
BOOTSTRAPCONFIDENCE INTERVAL CALCULATIONS
Based on 500bootstrap replicates
CALL :
boot.ci(boot.out= result)
Intervals :
Level      Normal              Basic
95%  (10.08, 11.30 )  (10.06, 11.26 )
Level    Percentile            BCa
95%  (10.20, 11.40 )  (10.13, 11.35 )
Calculationsand Intervals on Original Scale
Some BCaintervals may be unstable
传统估计：
2.5 %  97.5 %
(Intercept)31.20069 35.74810
eruptions  10.10996 11.34932
差别有些⼤，我们来看看qq图：
正态性不是很好，shapiro检验告诉我们⼏乎不可能认为是正态的。
bootstrap检验方法Shapiro-Wilk normality test
data:  faithful$waiting
W = 0.9221,p-value = 1.016e-10
观察下图，我们也可以看到waiting不服从正态分布，那么他的95%的置信区间可以通过以下代码获得：
mean1<-function(data,indices){
d<-data[indices,]
fit<-mean(d$waiting)
return(fit)
}
results<-boot(data=faithful,statistic=mean1,R=1000)
boot.ci(results)
输出结果：
BOOTSTRAPCONFIDENCE INTERVAL CALCULATIONS
Based on 1000bootstrap replicates
CALL :
boot.ci(boot.out= results)
Intervals :
Level      Normal              Basic
95%  (69.28, 72.50 )  (69.31, 72.53 )
Level    Percentile            BCa
95%  (69.26, 72.49 )  (69.10, 72.44 )
Calculationsand Intervals on Original Scale
与传统⽅法⽐较
[1] 69.2741872.51994
可见估计的稳健性。
⼆、  对多个统计量使⽤bootstrap⽅法
我们考虑博⽂⼆中的例⼦，获取⼀个统计向量（四个回归系数）的95%置信区间。      ⾸先，创建⼀个返回回归系数向量的函数：
betas<-function(formula,data,indices){
d<-data[indices,]
fit<-lm(formula,data=d)
return(fit$coef)
}
然后⾃助抽样500次：
results<-boot(data=states,statistic=betas,R=500,formula=Murder~.)
print(results)
输出结果：

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。