I'm interested in bias corrected bootstrap confidence intervals for my model predictions. I've written some code to perform the bootstrap, show below:
library(tidyverse)
#> Warning: package 'ggplot2' was built under R version 4.0.2
#> Warning: package 'tibble' was built under R version 4.0.2
#> Warning: package 'tidyr' was built under R version 4.0.2
#> Warning: package 'dplyr' was built under R version 4.0.2
library(boot)
x = rnorm(25)
eta = -0.5 + 0.2*x
y = rpois(25, exp(eta))
d = tibble(x, y)
predictions<-function(data, ind){
model = glm(y ~ x, data = data[ind,])
newdata = tibble(x = seq(-2,2, 0.01))
eta = predict(model, newdata = newdata)
exp(eta)
}
bt = boot(d, predictions, 1000)
boot.ci(bt, type = 'bca')
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 1000 bootstrap replicates
#>
#> CALL :
#> boot.ci(boot.out = bt, type = "bca")
#>
#> Intervals :
#> Level BCa
#> 95% ( 0.654, 6.095 )
#> Calculations and Intervals on Original Scale
#> Some BCa intervals may be unstable
The bootstrapping works as anticipated, however the boot.ci
call only computes a single confidence interval.
My desired output would be a dataframe which has:
- One column for bootstrap predictions (essentially, the mean of the bootstraps)
- One column for the lower
100*(1-a)
confidence limit, and - One column for the upper
100*(1-a)
confidence limit.
How can I use boot
to obtain the desired output?