Statisticians can reuse their data to quantify the uncertainty of complex models
Statistics is the branch of applied mathematics that studies ways of drawing inferences from limited and imperfect data. We may want to know how a neuron in a rat’s brain responds when one of its whiskers gets tweaked, or how many rats live in Manhattan, or how high the water will get under the Brooklyn Bridge, or the typical course of daily temperatures in the city over the year. We have some data on all of these things, but we know that our data are incomplete, and experience tells us that repeating our experiments or observations, even taking great care to replicate the conditions, gives more or less different answers every time. It is foolish to treat any inference from only the data in hand as certain.
If all data sources were totally capricious, there’d be nothing to do beyond piously qualifying every conclusion with “but we could be wrong about this.” A mathematical science of statistics is possible because, although repeating an experiment gives different results, some types of results are more common than others; their relative frequencies are reasonably stable. We can thus model the data-generating mechanism through probability distributions and stochastic processes—random series with some indeterminacy about how the events might evolve over time, although some paths may be more likely than others. When and why we can use stochastic models are very deep questions, but ones for another time. But if we can use them in a problem, quantities such as these are represented as “parameters” of the stochastic models. In other words, they are functions of the underlying probability distribution. Parameters can be single numbers, such as the total rat population; vectors; or even whole curves, such as the expected time-course of temperature over the year. Statistical inference comes down to estimating those parameters, or testing hypotheses about them.
These estimates and other inferences are functions of the data values, which means that they inherit variability from the underlying stochastic process. If we “reran the tape” (as Stephen Jay Gould used to say) of an event that happened, we would get different data with a certain characteristic distribution, and applying a fixed procedure would yield different inferences, again with a certain distribution. Statisticians want to use this distribution to quantify the uncertainty of the inferences. For instance, by how much would our estimate of a parameter vary, typically, from one replication of the experiment to another—say, to be precise, what is the root-mean-square (the square root of the mean average of the squares) deviation of the estimate from its average value, or the standard error? Or we could ask, “What are all the parameter values that could have produced this data with at least some specified probability?” In other words, what are all the parameter values under which our data are not low-probability outliers? This gives us the confidence region for the parameter—rather than a point estimate, a promise that either the true parameter point lies in that region, or something very unlikely under any circumstances happened—or that our stochastic model is wrong.
To get standard errors or confidence intervals, we need to know the distribution of our estimates around the true parameters. These sampling distributions follow from the distribution of the data, because our estimates are functions of the data. Mathematically the problem is well defined, but actually computing anything is another story. Estimates are typically complicated functions of the data, and mathematically convenient distributions all may be poor approximations of the data source. Saying anything in closed form about the distribution of estimates can be simply hopeless. The two classical responses of statisticians have been to focus on tractable special cases, and to appeal to asymptotic analysis, a method that approximates the limits of functions.
If you’ve taken an elementary statistics course, you were probably drilled in the special cases. From one end of the possible set of solutions, we can limit the kinds of estimator we use to those with a simple mathematical form—say, mean averages and other linear functions of the data. From the other, we can assume that the probability distributions featured in the stochastic model take one of a few forms for which exact calculation is possible, either analytically or via tables of special functions. Most such distributions have origin myths: The Gaussian bell curve arises from averaging many independent variables of equal size (say, the many genes that contribute to height in humans); the Poisson distribution comes from counting how many of a large number of independent and individually improbable events have occurred (say, radium nuclei decaying in a given second), and so on. Squeezed from both ends, the sampling distribution of estimators and other functions of the data becomes exactly calculable in terms of the aforementioned special functions.
That these origin myths invoke various limits is no accident. The great results of probability theory—the laws of large numbers, the ergodic theorem, the central limit theorem and so on—describe limits in which all stochastic processes in broad classes of models display the same asymptotic behavior. The central limit theorem (CLT), for instance, says that if we average more and more independent random quantities with a common distribution, and if that common distribution is not too pathological, then the distribution of their means approaches a Gaussian. (The non-Gaussian parts of the distribution wash away under averaging, but the average of two Gaussians is another Gaussian.) Typically, as in the CLT, the limits involve taking more and more data from the source, so statisticians use the theorems to find the asymptotic, large-sample distributions of their estimates. We have been especially devoted to rewriting our estimates as averages of independent quantities, so that we can use the CLT to get Gaussian asymptotics. Refinements to such results would consider, say, the rate at which the error of the asymptotic Gaussian approximation shrinks as the sample sizes grow.
To illustrate the classical approach and the modern alternatives, I’ll introduce some data: The daily closing prices of the Standard and Poor’s 500 stock index from October 1, 1999, to October 20, 2009. (I use these data because they happen to be publicly available and familiar to many readers, not to impart any kind of financial advice.) Professional investors care more about changes in prices than their level, specifically the log returns, the log of the price today divided by the price yesterday. For this time period of 2,529 trading days, there are 2,528 such values (see Figure 1). The “efficient market hypothesis” from financial theory says the returns can’t be predicted from any public information, including their own past values. In fact, many financial models assume such series are sequences of independent, identically distributed (IID) Gaussian random variables. Fitting such a model yields the distribution function in the lower left graph of Figure 1.
An investor might want to know, for instance, how bad the returns could be. The lowest conceivable log return is negative infinity (with all the stocks in the index losing all value), but most investors worry less about an apocalyptic end of American capitalism than about large-but-still-typical losses—say, how bad are the smallest 1 percent of daily returns? Call this number q0.01; if we know it, we know that we will do better about 99 percent of the time, and we can see whether we can handle occasional losses of that magnitude. (There are about 250 trading days in a year, so we should expect two or three days at least that bad in a year.) From the fitted distribution, we can calculate that q0.01=–0.0326, or, undoing the logarithm, a 3.21 percent loss. How uncertain is this point estimate? The Gaussian assumption lets us calculate the asymptotic sampling distribution of q0.01, which turns out to be another Gaussian (see the lower right graph in Figure 1), implying a standard error of ±0.00104. The 95 percent confidence interval is (–0.0347, –0.0306): Either the real q0.01 is in that range, or our data set is one big fluke (at 1-in-20 odds), or the IID-Gaussian model is wrong.
From its origins in the 19th century through about the 1960s, statistics was split between developing general ideas about how to draw and evaluate statistical inferences, and working out the properties of inferential procedures in tractable special cases (like the one we just went through) or under asymptotic approximations. This yoked a very broad and abstract theory of inference to very narrow and concrete practical formulas, an uneasy combination often preserved in basic statistics classes.
The arrival of (comparatively) cheap and fast computers made it feasible for scientists and statisticians to record lots of data and to fit models to them. Sometimes the models were conventional ones, including the special-case assumptions, which often enough turned out to be detectably, and consequentially, wrong. At other times, scientists wanted more complicated or flexible models, some of which had been proposed long before but now moved from being theoretical curiosities to stuff that could run overnight. In principle, asymptotics might handle either kind of problem, but convergence to the limit could be unacceptably slow, especially for more complex models.
By the 1970s statistics faced the problem of quantifying the uncertainty of inferences without using either implausibly helpful assumptions or asymptotics; all of the solutions turned out to demand even more computation. Perhaps the most successful was a proposal by Stanford University statistician Bradley Efron, in a now-famous 1977 paper, to combine estimation with simulation. Over the last three decades, Efron’s “bootstrap” has spread into all areas of statistics, sprouting endless elaborations; here I’ll stick to its most basic forms.
Remember that the key to dealing with uncertainty in parameters is the sampling distribution of estimators. Knowing what distribution we’d get for our estimates on repeating the experiment would give us quantities, such as standard errors. Efron’s insight was that we can simulate replication. After all, we have already fitted a model to the data, which is a guess at the mechanism that generated the data. Running that mechanism generates simulated data that, by hypothesis, have nearly the same distribution as the real data. Feeding the simulated data through our estimator gives us one draw from the sampling distribution; repeating this many times yields the sampling distribution as a whole. Because the method gives itself its own uncertainty, Efron called this “bootstrapping”; unlike Baron von Münchhausen’s plan for getting himself out of a swamp by pulling himself out by his bootstraps, it works.
Let’s see how this works with the stock-index returns. Figure 2 shows the overall process: Fit a model to data, use the model to calculate the parameter, then get the sampling distribution by generating new, synthetic data from the model and repeating the estimation on the simulation output. The first time I recalculate q0.01 from a simulation, I get -0.0323. Replicated 100,000 times, I get a standard error of 0.00104, and a 95 percent confidence interval of (–0.0347, –0.0306), matching the theoretical calculations to three significant digits. This close agreement shows that I simulated properly! But the point of the bootstrap is that it doesn’t rely on the Gaussian assumption, just on our ability to simulate.
The bootstrap approximates the sampling distribution, with three sources of approximation error. First there’s simulation error, using finitely many replications to stand for the full sampling distribution. Clever simulation design can shrink this, but brute force—just using enough replications—can also make it arbitrarily small. Second, there’s statistical error: The sampling distribution of the bootstrap reestimates under our fitted model is not exactly the same as the sampling distribution of estimates under the true data-generating process. The sampling distribution changes with the parameters, and our initial fit is not completely accurate. But it often turns out that distribution of estimates around the truth is more nearly invariant than the distribution of estimates themselves, so subtracting the initial estimate from the bootstrapped values helps reduce the statistical error; there are many subtler tricks to the same end. The final source of error in bootstrapping is specification error: The data source doesn’t exactly follow our model at all. Simulating the model then never quite matches the actual sampling distribution.
Here Efron had a second brilliant idea, which is to address specification error by replacing simulation from the model with resampling from the data. After all, our initial collection of data gives us a lot of information about the relative probabilities of different values, and in certain senses this “empirical distribution” is actually the least prejudiced estimate possible of the underlying distribution—anything else imposes biases or preconceptions, which are possibly accurate but also potentially misleading. We could estimate q0.01 directly from the empirical distribution, without the mediation of the Gaussian model. Efron’s “nonparametric bootstrap” treats the original data set as a complete population and draws a new, simulated sample from it, picking each observation with equal probability (allowing repeated values) and then re-running the estimation (as shown in Figure 2).
This new method matters here because the Gaussian model is inaccurate; the true distribution is more sharply peaked around zero and has substantially more large-magnitude returns, in both directions, than the Gaussian (see the top graph in Figure 3). For the empirical distribution, q0.01=–0.0392. This may seem close to our previous point estimate of –0.0326, but it’s well beyond the confidence interval, and under the Gaussian model we should see values that negative only 0.25 percent of the time, not 1 percent of the time. Doing 100,000 non-parametric replicates—that is, resampling from the data and reestimating q0.01 that many times—gives a very non-Gaussian sampling distribution (as shown in the right graph of Figure 3), yielding a standard error of 0.00364 and a 95 percent confidence interval of (–0.0477, –0.0346).
Although this is more accurate than the Gaussian model, it’s still a really simple problem. Conceivably, some other nice distribution fits the returns better than the Gaussian, and it might even have analytical sampling formulas. The real strength of the bootstrap is that it lets us handle complicated models, and complicated questions, in exactly the same way as this simple case.
To continue with the financial example, a question of perennial interest is predicting the stock market. Figure 4 is a scatter plot of the log returns on successive days, the return for today being on the horizontal axis and that of tomorrow on the vertical. It’s mostly just a big blob, because the market is hard to predict, but I have drawn two lines through it: a straight one in blue, and a curved one in black. These lines try to predict the average return tomorrow as functions of today’s return; they’re called regression lines or regression curves. The straight line is the linear function that minimizes the mean-squared prediction error, or the sum of the squares of the errors made in solving every single equation (called the least squares method). Its slope is negative (–0.0822), indicating that days with below-average returns tend to be followed by ones with above-average returns and vice versa, perhaps because people try to buy cheap after the market falls (pushing it up) and sell dear when it rises (pulling it down). Linear regressions with Gaussian fluctuations around the prediction function are probably the best-understood of all statistical models—their oldest forms go back two centuries now—but they’re more venerable than accurate.
The black curve is a nonlinear estimate of the regression function, coming from a constrained optimization procedure called spline smoothing: Find the function that minimizes the prediction error, while capping the value of the average squared second derivative. As the constraint tightens, the optimal curve, the spline, straightens out, approaching the linear regression; as the constraint loosens, the spline wiggles to try to pass through each data point. (A spline was originally a flexible length of wood craftsmen used to draw smooth curves, fixing it to the points the curve had to go through and letting it flex to minimize elastic energy; stiffer splines yielded flatter curves, corresponding mathematically to tighter constraints.)
To actually get the spline, I need to pick the level of the constraint. Too small, and I get an erratic curve that memorizes the sample but won’t generalize to new data; but too much smoothing erases real and useful patterns. I set the constraint through cross-validation: Remove one point from the data, fit multiple curves with multiple values of the constraint to the other points, and then see which curve best predicts the left-out point. Repeating this for each point in turn shows how much curvature the spline needs in order to generalize properly. In this case, we can see that we end up selecting a moderate amount of wiggliness; like the linear model, the spline predicts reversion in the returns but suggests that it’s asymmetric—days of large negative returns being followed, on average, by bigger positive returns than the other way around. This might be because people are more apt to buy low than to sell high, but we should check that this is a real phenomenon before reading much into it.
There are three things we should note about spline smoothing. First, it’s much more flexible than just fitting a straight line to the data; splines can approximate a huge range of functions to an arbitrary tolerance, so they can discover complicated nonlinear relationships, such as asymmetry, without guessing in advance what to look for. Second, there was no hope of using a smoothing spline on substantial data sets before fast computers, although now the estimation, including cross-validation, takes less than a second on a laptop. Third, the estimated spline depends on the data in two ways: Once we decide how much smoothing to do, it tries to match the data within the constraint; but we also use the data to decide how much smoothing to do. Any quantification of uncertainty here should reckon with both effects.
There are multiple ways to use bootstrapping to get uncertainty estimates for the spline, depending on what we’re willing to assume about the system. Here I will be cautious and fall back on the safest and most straightforward procedure: Resample the points of the scatter plot (possibly getting multiple copies of the same point), and rerun the spline smoother on this new data set. Each replication will give a different amount of smoothing and ultimately a different curve. Figure 5 shows the individual curves from 800 bootstrap replicates, indicating the sampling distribution, together with 95 percent confidence limits for the curve as a whole. The overall negative slope and the asymmetry between positive and negative returns are still there, but we can also see that our estimated curve is much better pinned down for small-magnitude returns, where there are lots of data, than for large-magnitude returns, where there’s little information and small perturbations can have more effect.
Smoothing Things Out
Bootstrapping has been ramified tremendously since Efron’s original paper, and I have sketched only the crudest features. Nothing I’ve done here actually proves that it works, although I hope I’ve made that conclusion plausible. And indeed sometimes the bootstrap fails; it gives very poor answers, for instance, to questions about estimating the maximum (or minimum) of a distribution. Understanding the difference between that case and that of q0.01, for example, turns out to involve rather subtle math. Parameters are functions of the distribution generating the data, and estimates are functions of the data or of the empirical distribution. For the bootstrap to work, the empirical distribution has to converge rapidly on the true distribution, and the parameter must smoothly depend on the distribution, so that no outlier ends up unduly influencing the estimates. Making “influence” precise here turns out to mean taking derivatives in infinite-dimensional spaces of probability distribution functions, and the theory of the bootstrap is a delicate combination of functional analysis with probability theory. This sort of theory is essential to developing new bootstrap methods for new problems, such as ongoing work on resampling spatial data, or model-based bootstraps where the model grows in complexity with the data.
The bootstrap has earned its place in the statistician’s toolkit because, of all the ways of handling uncertainty in complex models, it is at once the most straightforward and the most flexible. It will not lose that place so long as the era of big data and fast calculation endures.
- Efron, B. 1979. Bootstrap methods: another look at the jackknife. Annals of Statistics 7:1–26.