adjusted bootstrap confidence intervals (BCa) with parametric bootstrap in boot package

Question

I am attempting to use boot.ci from R's boot package to calculate bias- and skew-corrected bootstrap confidence intervals from a parametric bootstrap. From my reading of the man pages and experimentation, I've concluded that I have to compute the jackknife estimates myself and feed them into boot.ci, but this isn't stated explicitly anywhere. I haven't been able to find other documentation, although to be fair I haven't looked at the original Davison and Hinkley book on which the code is based ...

If I naively run b1 <- boot(...,sim="parametric") and then boot.ci(b1), I get the error influence values cannot be found from a parametric bootstrap. This error occurs if and only if I specify type="all" or type="bca"; boot.ci(b1,type="bca") gives the same error. So does empinf(b1). The only way I can get things to work is to explicitly compute jackknife estimates (using empinf() with the data argument) and feed these into boot.ci.

Construct data:

set.seed(101)
d <- data.frame(x=1:20,y=runif(20))
m1 <- lm(y~x,data=d)

Bootstrap:

b1 <- boot(d$y,
           statistic=function(yb,...) {
             coef(update(m1,data=transform(d,y=yb)))
           },
           R=1000,
           ran.gen=function(d,m) {
             unlist(simulate(m))
           },
           mle=m1,
           sim="parametric")

Fine so far.

boot.ci(b1)
boot.ci(b1,type="bca")
empinf(b1)

all give the error described above.

This works:

L <- empinf(data=d$y,type="jack",
            stype="i",
            statistic=function(y,f) {
              coef(update(m1,data=d[f,]))
            })

boot.ci(b1,type="bca",L=L)

Does anyone know if this is the way I'm supposed to be doing it?

update: The original author of the boot package responded to an e-mail:

... you are correct that the issue is that you are doing a parametric bootstrap. The bca intervals implemented in boot are non-parametric intervals and this should have been stated explicitely somewhere. The formulae for parametric bca intervals are not the same and depend on derivatives of the least favourable family likelihood when there are nuisance parameters as in your case. (See pp 206-207 in Davison & Hinkley) empinf assumes that the statistic is in one of forms used for non-parametric bootstrapping (which you did in your example call to empinf) but your original call to boot (correctly) had the statistic in a different form appropriate for parametric resampling.

You can certainly do what you're doing but I am not sure of the theoretical properties of mixing parametric resampling with non-parametric interval estimation.

Robbie

I perhaps have an advantage in that I have the book and it has a whole chapter on regression applications. I cannot tell, however, whether it would be useful to post my results of using jack.after.boot since you say you are convinced that the jackknife needs to be done first. – Fluent

Answer 1

After looking at the boot.ci page I decided to use a boot-object constructed along the lines of an example in Ch 6 of Davison and Hinkley and see whether it generated the errors you observed. I do get a warning but no errors.:

require(boot) 
lmcoef <- function(data, i){
      d <- data[i, ]
      d.reg <- lm(y~x, d)
      c(coef(d.reg)) }
lmboot <- boot(d, lmcoef, R=999)
m1
boot.ci(lmboot, index=2)   # I am presuming that the interest is in the x-coefficient
#----------------------------------
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = lmboot, index = 2)

Intervals : 
Level      Normal              Basic         
95%   (-0.0210,  0.0261 )   (-0.0236,  0.0245 )  

Level     Percentile            BCa          
95%   (-0.0171,  0.0309 )   (-0.0189,  0.0278 )  
Calculations and Intervals on Original Scale
Warning message:
In boot.ci(lmboot, index = 2) :
  bootstrap variances needed for studentized intervals

Fluent

hmm. thanks ... regression was just my example (my real problem is a GLMM), but I will take a look at this and try to figure out the difference between what you've done (which works) and what I've done (which doesn't) – Robbie

oops. now I see. The key is that I'm using parametric bootstrapping, you're not. I think if you re-do this with type="parametric" and the other necessary changes (i.e. adding the ran.gen function, etc.) that you'll end up with something that looks a lot like my example above and doesn't provide BCa estimates without semi-manual jackknife calculations ... – Robbie

What value does the type="parametric" provide in that construction? A linear (or linear on some other transformed scale) coefficient is a parametric estimate. So After reviewing the material in the 2nd chapter of Davison and Hinkley, I think I am doing a parametric bootstrap in the code above. – Fluent

I agree that the parametric analysis does not provide any advantage in the example. I am doing a generalized linear mixed model where the parametric bootstrap is definitely different from the standard nonparametric bootstrap, or even than a stratified bootstrap. – Robbie