What is Bootstrapping in Statistics? A Deep Dive

Bootstrapping confidence intervals for regression

Linear regression, the cornerstone of statistical modeling, is used to show the relationship between one or more independent variables and a dependent variable. When working with regression models, it’s common to assess uncertainty around our estimates. One way to do this is by calculating linear regression confidence intervals, which provide a range of values within which the true parameter, such as the mean, is likely to fall. Here, we will create confidence intervals using bootstrapping, and we will also create confidence intervals using the normal approximation method, also called the Wald confidence interval.

In the context of linear regression, bootstrap resampling involves randomly sampling from the dataset to create multiple bootstrap samples. We then fit a regression model to each sample. Finally, we use the distribution of model coefficients across these samples to estimate confidence intervals. In the following code, we use the bootstraps() function from the tidymodels package to perform bootstrap resampling on our dataset. We then create histograms to show the range of values that are present in our bootstrap resamples for the intercept and weight coefficients.

boots <- bootstraps(fish_df, times = 250, apparent = TRUE)

fit_lm_on_bootstrap <- function(split) {
    lm(Height ~ sqrt_Weight, analysis(split))
}

boot_models <-
   boots %>% 
   dplyr::mutate(model = map(splits, fit_lm_on_bootstrap),
         coef_info = map(model, tidy))

boot_coefs <- 
   boot_models %>% 
   unnest(coef_info)

percentile_intervals <- int_pctl(boot_models, coef_info)

ggplot(boot_coefs, aes(estimate)) +
   geom_histogram(bins = 30) +
   facet_wrap( ~ term, scales = "free") +
   labs(title="fish", subtitle = "Height ~ Square Root of Weight") +
   labs(caption = "95% confidence interval parameter estimates, intercept + estimate") +
   geom_vline(aes(xintercept = .lower), data = percentile_intervals, col = '#528fad') +
   geom_vline(aes(xintercept = .upper), data = percentile_intervals, col = '#528fad')

Linear model bootstrapped coefficients. Image by Author.

When we graph our possible regression lines, we will clip the most extreme lines on either end in order to get a 95% interval. The int_pctl() function makes this easy.

boot_aug <-
   boot_models %>%
   sample_n(50) %>%
   mutate(augmented = map(model, augment)) %>%
   unnest(augmented)

boots <- bootstraps(fish_df, times = 250, apparent = TRUE)

fit_lm_on_bootstrap <- function(split) {
    lm(Height ~ sqrt_Weight, analysis(split))
}

boot_models <-
   boots %>%
   dplyr::mutate(model = map(splits, fit_lm_on_bootstrap),
         coef_info = map(model, tidy))

boot_coefs <-
   boot_models %>%
   unnest(coef_info)

percentile_intervals <- int_pctl(boot_models, coef_info)

boot_aug <-
   boot_models %>%
   sample_n(50) %>%
   mutate(augmented = map(model, augment)) %>%
   unnest(augmented)

As a final step, let’s now create confidence intervals that would be generated from the more theoretical, closed-form equations. We can find the parameters by calling summary() of our linear model object.

linear_model <- lm(Height ~ sqrt(Weight), fish_df)
summary(linear_model)

Call:
lm(formula = Height ~ sqrt(Weight), data = fish_df)

Residuals:
    Min      1Q  Median      3Q     Max
-7.7425 -1.3611 -0.5628  2.2606  5.0131

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.31506    0.43636   5.305 3.77e-07 ***
sqrt(Weight)  0.37569    0.02186  17.183  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.533 on 157 degrees of freedom
Multiple R-squared:  0.6529, Adjusted R-squared:  0.6506
F-statistic: 295.3 on 1 and 157 DF,  p-value: < 2.2e-16

Here, the estimated coefficient for the sqrt_Weight variable is 0.37569. The standard error for sqrt_Weight is 0.02186. The degrees of freedom are 157 (which is 159 observations in our dataset minus two estimated parameters). And the t-value for a 95% confidence interval that is two-tailed is 1.975189. We know this final number through the qt() function: qt(0.975, 157). In the below visualization, the blue lines are the bootstrapped linear regression lines, and the dotted red lines are the confidence interval generated by the common analytical formula.

# Wald confidence interval
(CI_upper <- 0.37569 + 1.975 * 0.02186)
(CI_lower <- 0.37569 - 1.975 * 0.02186)

(CI_intercept_upper <- 2.31506 + 1.975 * 0.43636)
(CI_intercept_lower <- 2.31506 - 1.975 * 0.43636)

ggplot(boot_aug, aes(sqrt_Weight, Height)) +
   geom_line(aes(y = .fitted, group = id), alpha = .3, col = '#528fad') +
   geom_point(alpha = 0.005) +
   labs(title="fish", subtitle = "Height ~ Weight") +
   labs(caption = "linear regression bootstrap resampling") +
   geom_abline(intercept = 1.453249, slope = 0.4188635, linetype = "dashed", color = "#e76254") +
   geom_abline(intercept = 3.176871, slope = 0.3325165, linetype = "dashed", color = "#e76254")

Bootstrapped linear regressions with confidence intervals. Image by Author.

Here we can see clearly that our process of bootstrapping resampling has introduced a bit more uncertainty in our estimate. This is because the intervals for linear regression coefficients were calculated using theoretical formulas based on assumptions about the distribution of errors and the properties of the estimator. These formulas rely on assumptions such as normality of errors and constant variance.

Bootstrap resampling, on the other hand, is distribution-free, meaning that it makes minimal assumptions about the underlying data distribution. Instead, bootstrap resampling directly estimates the sampling distribution of our statistic of interest by resampling from the observed data. As a result, bootstrap confidence intervals can be more robust and reliable when the assumptions of traditional methods are violated or when dealing with small sample sizes, such as in the case of our Fish Market dataset with 159 observations.

Parametric vs. Non-Parametric Bootstrapping

In parametric bootstrapping, assumptions are made about the underlying distribution of the data, and resamples are generated based on those assumptions. This method is useful when you have prior knowledge or strong assumptions about the data’s distribution. Think about a sampling dataset that might not have a normal distribution, but the idea beyond the sampling is better known, so you can create a distribution using parameters from the population. 

Non-parametric bootstrapping, on the other hand, makes no assumptions about the data’s distribution. It resamples directly from the observed data with replacement, making it particularly valuable when the true distribution is unknown or hard to define. In our example above, we used non-parametric bootstrapping. Both methods allow you to estimate statistics such as standard errors and confidence intervals. However, non-parametric bootstrapping offers more flexibility for real-world datasets, especially when dealing with small or complex samples, and is more commonly used in practice.

Bootstrapping in Forecasting

In time series forecasting, bootstrapping can be applied to resample historical data and generate future forecasts, providing a distribution of possible outcomes rather than a single point estimate. This helps model the range of potential future scenarios and creates confidence intervals for predictions. Bootstrapping also underpins ensemble methods like bagging in time series models, which can reduce overfitting and improve the overall accuracy of the forecast by combining multiple models. Our Forecasting in R will teach you how to resample in time forecasting, whether you use ARIMA forecasting or another method.

Final Thoughts

I hope you have come to appreciate bootstrapping, if you didn’t already. Bootstrapping, as we have seen, is a powerful tool in both statistics and machine learning, and it offers a interesting way to estimate the variability and confidence of statistical measures without requiring strong assumptions about the underlying data. 

Consider starting our Machine Learning Scientist in Python career track and you sure to become an expert with working with distribution types and complex datasets. Our Statistical Inference in R skill track is another great option with a focus in hypothesis testing, randomization, and measuring uncertainty.