I'm using a wild bootstrap to explore the confidence intervals of a nonlinear regression mixed-effects model (specifically one that was solved using nlmer
). The model is a tad persnickety and although the nominal fit solves with no warnings or errors, I find that wild bootstrap replicates often fail to converge. For example, running 100 replicates yields a failure rate of about 50%.
What information is being lost as a result? Is there a generally accepted way to handle bootstrapping of a model that often fails to converge with replicated input data? Do the discarded non-convergent replicates compromise the conclusions from the other replicates?
Searching online has so far yielded no hits that talk about this regime. What I've found is along the lines of "individual bootstrap replications may occasionally fail to converge", e.g.: http://www.inside-r.org/packages/cran/hysteresis/docs/summary.ellipsefit
FWIW here's my R code structure.
n.bootstrap.samplings <- 100
n.bootstrap.current <- 0
n.bootstrap.failed <- 0
nlmem.fitted.values.bootstrapper <- function(data, indices){
#Wild bootstrap with Rademacher distribution
print('Wild bootstrap resampling...')
data.plot$target.value <- data.plot$value + scaling*data.plot$nlmem.resid*sample(c(-1, 1),
nrow(data.plot),
replace=T, prob=c(0.5,0.5))
return.value <- tryCatch({
#Refit the nlmer model on the resampled data
bs.nlmem <- perform.nlmer(data.plot)
fitted(bs.nlmem)
},
error = function(condition) {
assign('n.bootstrap.failed', n.bootstrap.failed + 1, envir = .GlobalEnv)
message('Error performing bootstrap fit.')
message(condition)
return (rep(NA, nrow(data.plot)))
}
)
assign('n.bootstrap.current', n.bootstrap.current + 1, envir = .GlobalEnv)
print(paste('Completed bootstrap replicate number', n.bootstrap.current))
return(return.value)
}
nlmem.fitted.values.bootstrapped <- boot(data=data.plot, statistic=nlmem.fitted.values.bootstrapper, R=n.bootstrap.samplings)
print('Bootstrap failure rate:')
print(n.bootstrap.failed/n.bootstrap.samplings)