The bootstrap is large sample approximation. Numerous confidence intervals have been developed from it.

For my data, which consists of fitting Generalized Additive Models (GAMs), I computed large-sample CIs for my parameter of interest. I then followed this with simulating the bootstrapped deviance residuals and adding them back into the model.

My issue is that sometimes the bootstrap returns strange CIs. Sometimes the same value is returned for both the upper and lower limits of the 95% CI.

Here is such an example using a cubic basis function (bs = "cr" in mgcv) with dimension k = 100. I also have 100 data points.


boot(data = boot.data, statistic = boot.fun, R = 1000)

Bootstrap Statistics :
    original       bias    std. error
t1* 28.52002 0.0002239734  0.00048061

Based on 1000 bootstrap replicates

boot.ci(boot.out = res, type = "all")

Intervals : 
Level      Normal              Basic         
95%   (28.52, 28.52 )   (28.52, 28.52 )  

Level     Percentile            BCa          
95%   (28.52, 28.52 )   (28.52, 28.52 )  
Calculations and Intervals on Original Scale
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable

The graphical output looks fine to me, but perhaps I should increase the number of resamples to 10,000 (but it's time consuming without parallelization).

enter image description here

I'm sure such a problem with bootstrap CIs has been encountered before, but I am just not sure of fixes or alternatives.

If it's of any help, here is my R code for the bootstrap.

## Bootstrapping to find Bias(N*) and SE(N*) ##

# Bootstrap setup

res <- resid(HAC.cr) - mean(resid(HAC.cr))  # center the residuals
n <- length(res)
boot.data <- data.frame(d,  fit = fitted(HAC.cr), res = res)
boot.fun <- function(data, i) {

# Updated GAM fit using the residuals

boot.fit <- gam(boot.data$means + res[i] ~ s(specs, bs = "cr", k = k), optimizer = c("outer", "bfgs"), data = data)

# Make sure the original estimate also gets returned

 if (all(i == 1:n)) {
    inv.predict(HAC.cr, y = p*max(HAC.mat), x.name = "specs", lower = lower, upper = upper)[1L] 
} else {
inv.predict(boot.fit, y = p*max(HAC.mat), x.name = "specs", lower = lower, upper = upper)[1L]

# Run bootstrap simulation

res <- boot(boot.data, boot.fun, R = 1000)  # collect bootstrap samples
plot(res) # histogram and QQplot
boot.ci(res, type = "all")  # obtain bootstrap confidence intervals

Note in the bootstrap function, I have used the integer data type [1L] over the standard numeric type (which is R's default). My reasoning for doing this is because R would always throw errors whenever [1L] was excluded.

Any assistance is kindly appreciated.

  • 2
    $\begingroup$ All the data is close to 28.52. So the confidence interval should be very tight and the lower endpoint and upper endpoint would be 28.52 because of rounding to two decimal places Maybe the results would be sensible if you carried a few more significant digits. $\endgroup$ – Michael Chernick Sep 27 '17 at 20:27
  • $\begingroup$ @MichaelChernick I will try that. I thought that perhaps there was an error in specifying my wrapper function for the bootstrap or in using the integer data type [1L]. $\endgroup$ – compbiostats Sep 27 '17 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.