Bootstrapping residuals: Am I doing it right?

First of all: From what I understood, bootstrapping residuals works as follows:

1. Fit model to data
2. Calculate the residuals
3. Resample the residuals and add them to 1.
4. Fit model to new dataset from 3.
5. Repeat n times, but always add the resampled residuals to the fit from 1.

Is that correct so far?

What I want to do is something slightly different:

I want to estimate parameter and prediction uncertainty for an algorithm that estimates some environmental variable.

What I have is a error-free time-series (from a simulation) of that variable, x_true, to which I add some noise, x_noise, in order to generate a synthetic dataset x. I then try to find optimal parameters by fitting my algorithm with the sum of squares sum((x_estimate - x_true)^2) (! not x_estimate - x !) as an objective function. In order to see how my algorithm performs and to create samples of my parameters' distributions, I want to resample x_noise, add it to x_true, fit my model again, rinse and repeat. Is that a valid approach to assess parameter uncertainty? Can I interpret the fits to the bootstrapped datasets as prediction uncertainty, or do I have to follow the procedure I posted above?

/edit: I think I haven't really made clear what my model does. Think of it as essentially something like a de-noising method. It's not a predictive model, it's an algorithm that tries to extract the underlying signal of a noisy time-series of environmental data.

/edit^2: For the MATLAB-Users out there, I wrote down some quick & dirty linear regression example of what I mean.

This is what I believe "ordinary" bootstrapping of residuals is (please correct me if I'm wrong): http://pastebin.com/C0CJp3d1

This is what I want to do: http://pastebin.com/mbapsz4c

• It will be clearer if you show the code that you have done so far. Aug 15, 2013 at 23:01
• I actually haven't coded anything so far in terms of bootstrapping. The code for my model is pretty complex, I don't thank that would help. As an example, we can assume that the model is a smoothing procedure like a moving average, with the moving window as the only model parameter. I have a series of (synthetic) measurements over time and add an error (not necessarily homoskedastic and normally distributed) to that. I then want to estimate the moving window which comes closest to the underlying "true" I know and want to assess uncertainty by bootstrapping my synthetic error. Does that help? Aug 15, 2013 at 23:09
• Here's some very bad MATLAB-style pseudo code, maybe it helps understand what I'd like to do: pastebin.com/yTRahzr5 Aug 15, 2013 at 23:24
• Sorry Fred, I don't know Matlab.Please tag as Matlab to get inputs from users. Aug 15, 2013 at 23:27
• Oh my question really isn't limited to MATLAB (and that isn't really MATLAB code, it's just some pseudo-code based on MATLABs syntax for for-loops and comments that wouldn't work anyway). But I can tag it just in case. Aug 15, 2013 at 23:32

Here is the general (semi-parametric-bootstrap) algorithm in more detail:

$\text{B}$ = number of bootstraps

the model:
$y = x\beta + \epsilon$

let $\hat{\epsilon}$ be the residuals

1. Run the regression and obtain the estimator(s) $\hat\beta$ and residuals $\hat\epsilon$.
2. Resample the residuals with replacement and obtain the bootstrapped residual vector $\hat\epsilon_\text{B}$.
3. Obtain the bootstrapped dependent variable by multiplying the estimator(s) from (1) with the original regressors and adding the bootstrapped residual: $y_\text{B} = x\hat\beta + \hat\epsilon_\text{B}$.
4. Run the regression with the bootstrapped dependent variables and the original regressors, this gives the bootstrapped estimator, i.e. regress $y_B$ on $x$, this gives $\hat\beta_\text{B}$.
5. Repeat the procedure $\text{B}$-times by going back to (2).
• Might be worth adding what this gets you over just having run step 1 and stopped (i.e. OLS). Is it in order to get an estimate on the variance of the parameters 𝛽? And what is the final estimate for 𝛽, the mean of the parameters found in all the runs?
– Dan
Aug 10, 2020 at 0:06
• The purpose is to get a confidence interval, say 95% of the mean estimate of B-hat using the distributions of B-hat from bootstrapping. The commenter should've noted that the # of bootstrapped residuals is equal to the dimensions of the y-vector, for example if your y-vector is 1000 predictions then the residual bootstrap vector is 1000 (from resampling w/replacement from the entire set of residuals from the regression) Oct 28, 2021 at 23:36

To see how an algorithm performs in terms of predictive accuracy/mean squared error, you probably need the Efron-Gong "optimism" bootstrap. This is implemented for easy use in the R rms package. See its functions ols, validate.ols, calibrate.

I'm not sure that my understanding is correct. But here is my suggestion to modify your code ("ordinary bootstrapping of residuals", lines 28-34) into:

for i = 2:n_boot
x_res_boot = x_residuals( randi(n_data,n_data,1) );
x_boot = x_res_boot+ x_best_fit;
p_est(:, i) = polyfit( t, x_boot, 1 );
x_best_fit2 = polyval( p_est(:, i), t );
x_residuals = x_best_fit2 - x_boot;
x_best_fit=x_best_fit2;
end


The idea is that each time you are using residuals not from the first run, but from the previous bootstrap fit. As for me, all other seems to be valid.

This is revised version that has been checked in MATLAB. Two errors have been fixed.

• Oh, well that was new to me. bsxfun is somewhat complicated; here is a new version that uses your idea and should be a bit clearer. However, that produces somewhat strange results. This is the result of always resampling the residuals of the first best fit and add them to the same (my initial idea), and this is what happens if I resample the residuals of each iteration and add them to each new best fit. Any ideas? Aug 19, 2013 at 12:33
• Whoops, small mistake in line 25 (should be p_est(:, i) instead of p_est(:, 1)), but even when I fix that the parameter distributions still look wonky: click Aug 19, 2013 at 12:47
• The answer is fixed and checked in MATLAB. Now it goes well. Aug 19, 2013 at 13:48
• New residuals for every fit - that was my first understanding of residual bootstrap. But I must admit that different sources are using residuals of fit to original data for that purpose. Here is nice tutorial on bootstrap (econ.pdx.edu/faculty/KPL/readings/mackinnon06.pdf). Seems that my approach is wrong while your implementation is right. Should I delete given answer? Aug 20, 2013 at 13:34
• Thanks for the follow up. IMHO, leave the answer for other users with the same question. I found that the literature (at least that which is available to me) isn't always clear on that subject and can be quite confusing. Aug 20, 2013 at 17:16

Spatiotemporal application

I applied @Sweetbabyjesus's answer. Unfortunately my model is more complex than an OLS regression (it's a distributed-lag 'hhh4 framework' spatiotemporal model from surveillance R package) causing instability when making such data changes, resulting in non-convergence messages.

If this happens to you then a second-best solution is producing jackknifed estimates by fitting the model to your dataset with one data point missing and repeating for $$N$$ to get a distribution of model estimates; in a spatiotemporal model context I interpret this as subtracting 1 from the response variable $$y$$ e.g. disease case counts while avoiding negative counts, so it is a trimmed subtraction function. Since I cannot fit the entire model if I make a single space x time cell go missing i.e. NA, subtracting by 1 count was the only option.

The mean of jackknifed estimates was very close to the model point estimates, however the jackknifed sd estimates underestimated the model sd by 1-3 orders of magnitude! Generally speak jackknifing is a 'first order approximation' to a bootstrap, hence why I see this solution as suboptimal.