2
$\begingroup$

I'm interested in using the weighted bootstrap to correct for selection bias with a known form. I simulated a very simple example where the underlying data, $X$, are $N(0,1)$ and we are calculating a sample mean. However, there is selection bias such that every $X>0$ gets sampled, while only 5% of $X<0$ get sampled. Thus, the sample mean is biased.

To correct for this, I resampled with replacement using inverse-probability-of-selection weights. I did this using the package boot with weights equal to 1 or 20. I computed confidence intervals using either the percentile method or the raw empirical quantiles of the bootstraps. Even though the resampled means were unbiased, both types of confidence interval had awful coverage (~53%).

As a sanity check, I set the parameter that controls the ratio of selection probabilities (eta below) to 1, such that there is no selection bias and we are just doing plain-vanilla resampling. Lo and behold, the coverage then is fine (96.5%). Since the bootstrap works for these data without selection, I don't think the excellent answers here hold (e.g., issues with non-pivotality, asymptotics, etc.).

Also, although I realize other methods of CI construction (e.g., BCa or basic) are sometimes better, I can't use them in this case because the original sample estimate is obviously biased so can't be used as a benchmark.

My code is below:

##### Helper Fn: Extract CI Limits from boot.ci #####
# list with first entry for b and second entry for t2
# n.ests: how many parameters were estimated?
get_boot_CIs = function(boot.res, type, n.ests) {
  bootCIs = lapply( 1:n.ests, function(x) boot.ci(boot.res, type = type, index = x) )

  # list with first entry for b and second entry for t2
  # the middle index "4" on the bootCIs accesses the stats vector
  # the final index chooses the CI lower (4) or upper (5) bound
  bootCIs = lapply( 1:n.ests, function(x) c( bootCIs[[x]][[4]][4],
                                             bootCIs[[x]][[4]][5] ) )
}


##### Helper Fn: Check CI coverage #####
covers = function( truth, lo, hi ) {
  return( (lo <= truth) & (hi >= truth) )
}


sim.reps = 200
boot.reps = 500
n = 5000  # initial sample size before selection
library(boot)

# sanity check: using eta = 1 gives 96.5% coverage for both methods


for ( i in 1:sim.reps ) {

    x = rnorm(n=5000)

    # sample values non-randomly
    # positive ones have an eta-fold higher chance of being sampled
    eta = 20
    weight = rep(1, length(x))
    weight[x < 0] = eta

    # indicator for whether we sample each observation
    keep = rbinom(n=n, size=1, prob=1/weight)

    d = data.frame(x, weight)
    d = d[ keep == 1, ]

    # this will be > 0
    # mean(d$x)

      boot.res = boot( data = d, 
                       parallel = "multicore",
                       R = boot.reps, 
                       weights = d$weight,
                       statistic = function(original, indices) {

                         b = original[indices,]
                         mean(b$x)
                       }
                      )  # end call to boot()

      # in case there is some weird problem with
      #  computing the boot CI
      tryCatch( {
        percCIs = get_boot_CIs(boot.res, "perc", n.ests = 1)
      }, error = function(err) {
        percCIs <<- list( c(NA, NA), c(NA, NA) )
      } )

      # quantiles
      qlo = quantile(boot.res$t[,1], 0.025)
  qhi = quantile(boot.res$t[,1], 0.975)

    rows =     data.frame( 
                            Method = c( "PercBT",
                                        "JustQuantiles" ),

                            # stats for mean estimate
                            # note that both boot CI methods use the same point estimate
                            XBarWtd = c( mean( boot.res$t[,1] ),
                                   mean( boot.res$t[,1] ) ),

                            Lo = c( percCIs[[1]][1],
                                    qlo ),

                            Hi = c( percCIs[[1]][2],
                                    qhi ),

                            Coverage = c( covers( 0, percCIs[[1]][1], percCIs[[1]][2] ),
                                          covers( 0,
                                                  qlo,
                                                  qhi ) )
    )

    if ( i == 1 ) res = rows
    else res = rbind(res, rows)
}

# see the underwhelming results
library(dplyr)
vars = c("XBarWtd", "Lo", "Hi", "Coverage")
res %>% group_by(Method) %>% summarise_at( vars(vars), mean )

And the disappointing results:

  Method        XBarWtd      Lo     Hi Coverage
  <fct>           <dbl>   <dbl>  <dbl>    <dbl>
1 JustQuantiles 0.00430 -0.0333 0.0423    0.53 
2 PercBT        0.00430 -0.0341 0.0430    0.535
$\endgroup$
3
$\begingroup$

The weights argument in the boot function is looking for resampling importance weights, not inverse probability weights. When no weights are specified, the bootstrap assumes an importance resampling weight of $n^{-1}$ for each unit, where $n$ is the sample size (i.e., uniform weighting). Importance weights are resampling probabilities that are intended to improve the efficiency of Monte Carlo estimates of bootstrap quantiles. See, e.g., this reference: "On importance resampling for the bootstrap" for a description of importance weights and their intended purpose. As to your problem, if you look at the boot function's source code, you will see the following:

if (!is.null(weights)) 
  weights <- t(apply(matrix(weights, n, length(R), 
                            byrow = TRUE), 2L, normalize, strata))

And digging a little deeper, you find that

boot:::normalize

normalizes the weights to add to 1. So, essentially, when you provide a vector of 1's and 20's, it becomes a vector of 1/sum(vector)'s and 20/sum(vector)'s. Boot then oversamples the 20/sum(vector)'s 20 times more than the 1/sum(vector)'s, so it is a happy coincidence that this is precisely the (weighted) correction needed to make your estimate of the mean unbiased. It is not, however, the intention of the weights argument, so it makes sense that other estimated quantities are not correct.

$\endgroup$
  • $\begingroup$ Thanks; good point. However, I tried writing my own weighted resampling (using dplyr::sample_n) and calculating either the raw quantiles or the basic CI myself (the latter using an unbiased weighted average from the original data), and both methods still had very poor coverage. Any thoughts? $\endgroup$ – half-pass Jan 2 '19 at 13:29
  • $\begingroup$ Seems like I answered your original question, which was something like "why is CI coverage bad when I use Horvitz-Thompson weights with package boot?" Now you have a new question, which is something like, "how does one account for probability sampling weights in a nonparametric bootstrap?" That's an interesting question that deserves a new post with some code showing what you have done. And also, no upvote? I took the time to read your code, replicate your problem, dig into the boot package and write up an answer with a relevant reference. I'm feeling a little burned over here. $\endgroup$ – bsbk Jan 3 '19 at 0:20
  • $\begingroup$ Okay, I apologize for not upvoting earlier. However, what I was trying to say is that I believe the explanation you provided -- while a very good guess -- is not actually the reason for the bad behavior, because the experiment I describe in the comment appears to falsify the explanation. Does that make sense? $\endgroup$ – half-pass Jan 3 '19 at 0:45
  • $\begingroup$ In your original post, above, you expected correct CI coverage because you were passing probability weights to boot. In my answer, I detailed why that was an unreasonable expectation and, in doing so, answered the question. Your new question is about reasonable approaches for variance estimation when bootstrapping with sampling weights. It is deserving of a new post in which you give the code you wrote that is performing poorly and related references (if any) that motivate your approach. I suggest you look at Lumley's Survey package before posting. $\endgroup$ – bsbk Jan 3 '19 at 16:20

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.