{"id":213,"date":"2018-06-13T21:19:58","date_gmt":"2018-06-13T21:19:58","guid":{"rendered":"http:\/\/statisticalbiophysicsblog.org\/?p=213"},"modified":"2018-12-09T05:08:20","modified_gmt":"2018-12-09T05:08:20","slug":"recovering-from-bootstrap-intoxication","status":"publish","type":"post","link":"https:\/\/statisticalbiophysicsblog.org\/?p=213","title":{"rendered":"Recovering from bootstrap intoxication"},"content":{"rendered":"

\n I want to talk again<\/a> today about the essential topic of analyzing statistical uncertainty – i.e., making error bars – but I want to frame the discussion in terms of a larger theme: our community\u2019s often insufficiently critical adoption of elegant and sophisticated ideas. I discussed this issue a bit previously in the context of PMF calculations<\/a>. To save you the trouble of reading on, the technical problem to be addressed is statistical uncertainty for high-variance data with small(ish) sample sizes.\n<\/p>\n

\n \n<\/p>\n

\n Elegance and effectiveness are not the same thing. I should know something about it because I used to try to make a living in part by pursuing algorithmic elegance (\u2018resolution exchange\u2019 simulation<\/a> and path sampling for Jarzynski-based free energy calculations<\/a>, to name two of my group\u2019s less-than-world-changing developments). In other words, I\u2019m both a perpetrator and victim of the problem.\n<\/p>\n

\n Let\u2019s dive into the bootstrap statistical analysis approach, which is truly beautiful and mathematically well-founded, but which can fail in at least one very important way – if we aren\u2019t careful about using it in a manner consistent with its underlying assumptions<\/em>. That\u2019s an academic way of saying that I was trying to use bootstrapping \u2018off the shelf\u2019 without understanding it well. And in some sense, the problem I will discuss results from insufficient data even though the bootstrap is sometimes noted as useful for small samples \u2026 but also as potentially unreliable in that regime<\/a>; see Schenker as well as Chernick & Labudde articles listed below.\n<\/p>\n

\n First, what is the basic bootstrapping idea? It is a perfectly valid strategy for estimating statistical variation of arbitrary observables based on a given sample of data. The sample at hand (e.g., a set of configurations, or a set of energy values or rates) is assumed to be representative and hence is used as a proxy for the true distribution. From this proxy distribution, which is discrete even if the underlying space is continuous, numerous \u201cbootstrap\u201d resamples<\/a> can be drawn, which are just samples with replacement. These bootstrap (re)samples represent a simulation of all possible samples assuming<\/em> the initial sample is representative – and hence represent all variation consistent with that assumption.\n<\/p>\n

\n With a set of bootstrap samples in hand, it is straightforward to make a confidence interval for an arbitrary <\/em>observable. For each sample we have, we can calculate that observable and make a histogram. For example, if each sample consists of a single scalar \"x\", we could compute the \u201cobservable\u201d \"x^3\" for every sample and build a histogram of these. Likewise, from a set of configurations, we could build a histogram of energy values. Then, a 95% confidence interval would consist of the 2.5 and 97.5 %ile values from the histogram.\n<\/p>\n

\n For later, remember this bootstrap confidence interval represents the variation expected if the original sample was a good proxy for the true distribution.<\/em> This is a \u201cfrequentist\u201d idea that we\u2019ll revisit below. For now it sounds reasonable enough. And of course, bootstrapping is a very elegant procedure!\n<\/p>\n

\n Let\u2019s move on to the downside. We\u2019ll consider a concrete, but very simple case where a failure can occur. Before you read on, or tell yourself a \u201ctoy\u201d example is not important, I want to emphasize that our simple binomial example is highly relevant to certain types of real data, as I will describe later on in the post.<\/em>\n<\/p>\n

\n Assume we want to estimate the probability \"p\" of the value 1 in a binomial distribution (where zero occurs with probability \"1-p\"). That is, \"p\" itself is the observable of interest. Say the sample we have in hand to estimate \"p\" is (1,0,0,0,0). Of course we would estimate \"p=1/5\". But how confident are we based on just five draws from the binomial distribution? Perhaps \"p\" is actually much smaller and we were lucky to see a single 1. Or perhaps \"p\" was larger and we got unlucky. How do we think about this?\n<\/p>\n

\n The bootstrapping procedure readily provides a confidence interval, as described above \u2026 but we will see that it is seriously flawed. Bootstrapping calls for sampling with replacement repeatedly from the original (1,0,0,0,0) sample. Sampling with replacement means that we draw sets of five elements from the original set and every element is chosen with equal probability. Thus, any element may randomly occur more than once (which is probably easier to think about if we had five different values, but we\u2019ll stick to binomial). In our binomial case, a couple of the (re)samples might look like (0,1,0,0,1) or (0,0,0,0,0). Imagine drawing thousands of such samples and estimating \"p\" from each, simply based on the fraction of ones \u2013 i.e., 2\/5 or 0 for the preceding examples. Then we can construct a histogram of these values and use the 2.5 and 97.5%iles to make a confidence interval. This interval represents the range of the 95% most likely values if<\/em> our original sample was representative of the true distribution.\n<\/p>\n

\n For the binomial case, there are few enough possibilities that we can easily predict what will happen in a bootstrap process. In particular, we know that in generating each bootstrap sample we have a 1\/5 chance of a 1 and 4\/5 chance of 0. Let\u2019s focus on a key occurrence, the outcome (0,0,0,0,0), which leads to \"p=0\" and has a probability of \"(4/5)^5 \simeq 0.33\". This tells us right away that the lower limit of the confidence interval will be 0 (because the 2.5% is less than 33% and hence occurs at \"p=0\" in the histogram). In fact, since \"p=0\" occupies 33% of the probability, even when we chop off the lowest 2.5%, bootstrapping\u2019s confidence interval still covers for 30% of its 95%!\n<\/p>\n

\n And there we have the problem: the bootstrapping confidence interval wrongly indicates that there\u2019s a decent chance that \"p=0\". Why is that wrong? Well, our sample of just five binomial values already has a 1 in it. Although we don\u2019t have much data, we do have enough data to rule out certain models:<\/em> \"p\" must be finite and not zero. We can use Bayesian ideas to quantify our common sense, but let\u2019s hold off on that for a moment.\n<\/p>\n

\n Let me now explain why a binomial example has been relevant for my own research on reasonably complex biomolecular systems. We have been simulating the folding of proteins using the \u201cweighted ensemble\u201d (WE) path sampling method, which can yield unbiased rate estimates. Each (quite expensive) WE run yields an estimate for the folding rate, but there is a very large variance among these estimates: one run may yield a rate that is several orders of magnitude larger than another. And a set of 10-30 WE folding runs yields a set of rates whose average is dominated by a small number (< 5) of the runs. The small number of large values is essentially like the 1\u2019s in a binomial distribution, with the remainder like the 0\u2019s. In that sense it\u2019s the same situation, and bootstrapping will suggest an artifactually small lower limit for the confidence interval of the folding rate. So for me, the binomial example is very relevant. (Note that a \u201cregular\u201d confidence interval based on the standard error of the mean will yield unphysical negative values because the data is highly non-Gaussian. Note also that one difference between the folding data and binomial samples is that in the binomial case we know that 1 is the largest possible value whereas we don\u2019t know the largest possible value for the folding rate ahead of time. However, the critique of the lower confidence limit is still valid.)\n<\/p>\n

\n We\u2019ve seen what goes wrong procedurally with the bootstrap in for a high-variance small sample, but what goes wrong theoretically? Here we start to get into issues of frequentism and Bayesianism, but I want to avoid jargon and abstraction. Let\u2019s just remind ourselves what the bootstrap procedure does. It generates a rather complete distribution of outcomes assuming that the sample in hand is a good approximation to the true underlying distribution.<\/em> But in generating a confidence interval for, say, the \"p\" parameter of a binomial distribution, I would argue it\u2019s clear that it\u2019s not really the possible outcomes we care about. Rather, we want to estimate the likelihood that a given \"p\" could have generated the sample we have, which FYI is the Bayesian point of view: find the models most consistent with data. Also, as with many statistical approaches, bootstrapping seems to be theoretically founded on large-numbers statistics, which clearly are lacking in our example.\n<\/p>\n

\n More concretely, given the binomial sample (1,0,0,0,0), it\u2019s clear this is consistent with \"p=1/5\", but what about other values of \"p\"? This is easy to quantify in the binomial case because from elementary binomial statistics<\/a>, the probability of a sample with a single 1 is \"p_{\mathrm{sample}} = 5p(1-p)^4\". From a Bayesian perspective, we can use this expression as a basis for estimating the probability that our sample came from any<\/em> \"p\". That is, even though strictly speaking the expression for \"p_{\mathrm{sample}}\" gives the fraction of five-element samples expected to have a single 1 for a fixed value<\/em> of \"p\", we can reverse the direction of the logic and use \"p_{\mathrm{sample}}\" as our best guess for the probability (density) that our particular sample came from any given value of \"p\". So aside from the issue of normalization, we can think of it as a function of \"p\": \"p_{\mathrm{sample}}(p)\" . Notice that this function has a maximum at \"p=1/5\" and it correctly tells us there is zero<\/em> probability that \"p=0\" or \"p=1\", in sharp contrast to bootstrapping.\n<\/p>\n

\n That\u2019s where we have to end this discussion. I don\u2019t have the time, space, or expertise to give a complete description of the best way to proceed with small-sample data of a binomial character. (I think a partial solution involves \u2018Bayesian bootstrapping\u2019 – see our initial efforts here<\/a>.) But I hope you have understood the basic reasons why even a beautiful idea like bootstrapping can go wrong. Of course, it\u2019s fair to say that I was guilty of trying to use bootstrapping without understanding it thoroughly enough. In any case, I hope my experience will be useful for you.\n<\/p>\n

\n In short, question authority, question precedent! Examine the assumptions underlying the calculations you do.\n<\/p>\n

Further reading<\/strong>\n<\/p>\n