How to deal with bootstrap replicates that fail to converge?

Question

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)

Answer 1

You are going to wind up with anti-conservative estimates of the errors of parameters in the non-linear mixed effects model. In larger samples, we presume that your model would behave more stably. These bootstrap replicates which are diverging to $\infty$ would necessarily be larger estimates than those in the normal range of estimates which are properly converging. To convince yourself of this, you can take the simple issue of separability in bootstrap estimation of confidence regions for logistic regression. See the example code below which honestly has kinks, but can be modified to answer this question.

Bootstraps make some assumptions about the distribution of your test statistic. If there is a non-zero probability weight at $\infty$ according to the bootstrapped distribution of parameter estimates, you should consider a better way of creating resampling based statistics, such as a parametric bootstrap, or possibly a permutation test.

set.seed(1)
n <- 8
x <- seq(n)

do.one <- function(n) {
  x <- seq(n)
  y <- rbinom(n, 1, x/n)

  b <- glm.fit(cbind(1,x), y, family=binomial())$coef

  boot.sd <- replicate(100, {
    x <- sample(x, replace=T)
    glm.fit(cbind(1,x), y[x], family=binomial())$coef
  })

  ## diverged if coef greater than
  max.div <- 10000
  keep <- apply(abs(boot.sd), 2, max) < max.div
  boot.sd <- boot.sd[, keep]

  c(b, apply(boot.sd, 1, sd))
}

options(warn=-1) ## many warning on divergent values
a <- replicate(100, do.one(10))

## actual sampling distribution of mean:
apply(a[1:2, ], 1, sd)

## estimated sampling distribution of mean
apply(a[3:4, ], 1, mean)