经管之家App
让优质教育人人可得
立即打开
Bootstrap is one of the most famous resampling technique and is very useful to get confidence intervals in situations where classical approach (t- or z- tests) would fail.
What is bootstrapInstead of writing down some equations let’s directly see how one may perform bootstrap. Below we will show a simple bootstrap example using the height of 100 person in the population.
set.seed(20151101)height<-rnorm(100,175,6)Now we will resample with replacement 1000 times and compute the median:
t0<-median(height)t<-sapply(1:1000,function(x) median(sample(x=height,size=100,replace=TRUE)))hist(t)abline(v=t0,col="orange",lwd=3)Here is the histogram we get:
And that’s it, this is the essence of bootstrap: resampling the observed data with replacement and computing the statistic of interest (here the median) many times on the resampled data to get a distribution of the statistic of interest. This distribution of the statistic of interest can then be used to compute, for example, confidence intervals.
When to use bootstrapBootstrap is used to enable inference on the statistic of interest when the true distribution of this statistic is unknown. For example in linear model the parameter of interest have a known distribution from which standard errors and formal tests can be performed. On the other hand for some statistics (median, differences between two models …), if the analyst do not want to spend time writing down equations,bootstrapping might be a great approach to get standard errors and confidence intervals from the bootstrapped distribution.
When will bootstrap failThere are some situations where bootstrapped will fail: (i) the statistic of interest is at the edge of the parameter space (like minimum or maximum), the bootstraped distribution does not converge (as the number of bootstrap sample increase to infinity) to the true distribution of the statistic. (ii) sample size is small, bootstrapping will not increase the power of statistical tests. If you sample to few data to detect an effect of interest, using bootstrap will not magically solve your problem even worse the bootstrap approach will perform less well than others.
How many bootstrap samplesAs much as possible will be the answer. Note that in this post I used low numbers to speed up the computations on my “old” computer.
Non-parametric and parametric bootstrap using the boot libraryThe boot library in R is very convenient to easily compute confidence intervals from bootstrap samples. Non-parametric bootstrap with boot:
library(boot)b1<-boot(data=height,statistic=function(x,i) median(x),R=1000)boot.ci(b1)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicatesCALL : boot.ci(boot.out = b1)Intervals : Level Normal Basic 95% (174.1, 176.2 ) (174.2, 176.4 ) Level Percentile BCa 95% (173.9, 176.1 ) (173.8, 176.1 ) Calculations and Intervals on Original ScaleParametric bootstrap with boot:
x<-runif(100,-2,2)y<-rnorm(100,1+2*x,1)dat<-data.frame(x=x,y=y)Simple example with a linear model
m<-lm(y~x)We are interested in getting the confidence intervals for the coefficient of the model:
foo<-function(out){ m<-lm(y~x,out) coef(m)}The function rgen generate new response vector from the model:
rgen<-function(dat,mle){ out<-dat out$y<-unlist(simulate(mle)) return(out)}Generate 1000 bootstrap sample
b2<-boot(dat,foo,R=1000,sim="parametric",ran.gen=rgen,mle=m)Compute the confidence intervals for the two coefficients:
boot.ci(b2,type="perc",index=1)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicates CALL : boot.ci(boot.out = b2, type = "perc", index = 1)Intervals : Level Percentile 95% ( 0.8056, 1.1983 ) Calculations and Intervals on Original ScaleCompute the confidence intervals for the two coefficients:
boot.ci(b2,type="perc",index=2)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 1000 bootstrap replicatesCALL : boot.ci(boot.out = b2, type = "perc", index = 2)Intervals : Level Percentile 95% ( 1.871, 2.217 ) Calculations and Intervals on Original ScaleIn the non-parametric case boot expect two arguments in the function returning the statistic of interest: the first one is the object from which to compute the statistic and the second is a vector of index (i), frequencies (f) or weight (w) defining the bootstrap sample. In our example, boot will generate a series of indices (named i) with replacement and these will be used to subset the original height vector. In the parametric case the function returning the statistic(s) of interest only need one argument: the original dataset. We then need to supply another function (ran.gen) describing how to generate the new data, it needs to return an object of the same form as the original dataset. This random data generating function need two arguments: the first one is the original dataset and the second one contain maximum likelihood estimate for the parameter of interest, basically a model object. The new dataset generated by theran.gen function will then be passed to the statistic function to compute the bootstrapped value for the statistic of interest. It is then straightforward to get the confidence intervals for the statistic usingboot.ci.
Bootstrap applied to mixed-effect modelsMixed-effect models are rather complex and the distributions or numbers of degrees of freedom of various output from them (like parameters …) is not known analytically. Which is why the author of the lme4 package recommend the use of bootstrap to get confidence intervals around the model parameters, the predicted values but also to get p-values from likelihood ratio tests.
A simple random intercept model:
dat<-data.frame(x=runif(100,-2,2),ind=gl(n=10,k=10))dat$y<-1+2*dat$x+rnorm(10,0,1.2)[dat$ind]+rnorm(100,0,0.5)m<-lmer(y~x+(1|ind),dat)Get the bootstrapped confidence intervals for the model parameters:
b_par<-bootMer(x=m,FUN=fixef,nsim=200)boot.ci(b_par,type="perc",index=1)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 200 bootstrap replicatesCALL : boot.ci(boot.out = b_par, type = "perc", index = 1)Intervals : Level Percentile 95% ( 0.393, 1.824 ) Calculations and Intervals on Original ScaleSome percentile intervals may be unstableGet the bootstrapped confidence intervals for the model parameters:
boot.ci(b_par,type="perc",index=2)BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONSBased on 200 bootstrap replicatesCALL : boot.ci(boot.out = b_par, type = "perc", index = 2)Intervals : Level Percentile 95% ( 1.857, 2.041 ) Calculations and Intervals on Original ScaleSome percentile intervals may be unstableOr alternatively:
confint(m,parm=c(3,4),method="boot",nsim=200,boot.type="perc") 2.5 % 97.5 %(Intercept) 0.4337793 1.819856x 1.8611089 2.058930Get the bootstrapped confidence intervals around the regression curves:
new_dat<-data.frame(x=seq(-2,2,length=20))mm<-model.matrix(~x,data=new_dat)predFun<-function(.) mm%*%fixef(.)bb<-bootMer(m,FUN=predFun,nsim=200) #do this 200 timesAs we did this 200 times the 95% CI will be bordered by the 5th and 195th value:
bb_se<-apply(bb$t,2,function(x) x[order(x)])new_dat$LC<-bb_se[1,]new_dat$UC<-bb_se[2,] new_dat$pred<-predict(m,newdata=new_dat,re.form=~0)Plot the results
plot(y~x,dat,pch=16)lines(pred~x,new_dat,lwd=2,col="orange")lines(LC~x,new_dat,lty=2,col="orange")lines(UC~x,new_dat,lty=2,col="orange")Here is the plot:
Finally get bootstrapped p-values from the likelihood ratio test between two models.
library(pbkrtest)m_0<-update(m,.~.-x)PBmodcomp(m,m_0,nsim=200)Parametric bootstrap test; time: 14.35 sec; samples: 200 extremes: 0;large : y ~ x + (1 | ind)small : y ~ (1 | ind) stat df p.value LRT 275.19 1 < 2.2e-16 ***PBtest 275.19 0.004975 ** ---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Drawing confidence intervals around the regression curves is tricky due to the random effect estimated values which comes with there own uncertainty (have a look at dotplot(ranef(m,condVar=TRUE)) to see it). Bootstrapping is an efficient way to take these uncertainties into account since the random deviates are re-computed for each draw. Finally getting p-values for the effect of a fixed-effect term can be done using a parametric bootstrap approach as described here and implemented in the function PBmodcompfrom the pbkrtest package. In the output of PBmodcomp the bootstrapped p-values is in the PBtest line, the LRT line report the standard p-value assuming a chi-square distribution for the LRT value.
ConclusionWith computers being always faster, bootstrap enable us to get reliable confidence interval estimate (given that your original sample size is large enough) without relying on some hypothetical distributions. Bootstrap can also be extended in a bayesian framework, see for example this nice post.
If you have any comments or suggestion feel free to post a comment below.
扫码加好友，拉您进群
全部版块
我的主页
收藏
