6.10 Bayesian bootstrap regression
We implement the Bayesian bootstrap (Donnald B. Rubin 1981) for linear regression models. In particular, the Bayesian bootstrap simulates the posterior distributions by assuming that the sample cumulative distribution function (CDF) is the population CDF (this assumption is also implicit in the frequentist bootstrap (B. Efron 1979)).
Given \(y_i \stackrel{i.i.d.}{\sim} \mathcal{F}\), where \(\mathcal{F}\) does not specify a particular parametric family of distributions, but instead sets \(\mathbb{E}(y_i \mid \boldsymbol{x}_i) = \boldsymbol{x}_i^{\top} \boldsymbol{\beta}\), with \(\boldsymbol{x}_i\) being a \(K\)-dimensional vector of regressors and \(\boldsymbol{\beta}\) a \(K\)-dimensional vector of parameters, the Bayesian bootstrap generates posterior probabilities for each \(y_i\), where the values of \(\boldsymbol{y}\) that are not observed have zero posterior probability.
The algorithm to implement the Bayesian bootstrap is the following:
Algorithm: Bayesian Bootstrap
Draw g from a Dirichlet distribution
\[ \mathbf{g} \sim \text{Dir}(\alpha_1, \alpha_2, \ldots, \alpha_N), \] with \(\alpha_i = 1\) for all \(i\).The vector
\[ \mathbf{g} = (g_1, g_2, \ldots, g_N) \] contains the probability weights assigned to each observation
\((y_1,\mathbf{x}_1), (y_2,\mathbf{x}_2), \ldots, (y_N,\mathbf{x}_N)\)
for each Bayesian bootstrap replication.For \(s = 1, \ldots, S\):
- Sample the pairs \((y_i,\mathbf{x}_i)\) N times with replacement using probabilities \(g_i\), \(i = 1,\ldots,N\).
- Estimate \(\boldsymbol{\beta}^{(s)}\) using weighted least squares in the model
\[ E(\mathbf{y}\mid \mathbf{X}) = \mathbf{X}\boldsymbol{\beta}, \] where the regression weights are the Dirichlet weights \(g_i\).
The posterior distribution of \(\boldsymbol{\beta}\) under the Bayesian bootstrap is the empirical distribution of
\[ \boldsymbol{\beta}^{(1)}, \boldsymbol{\beta}^{(2)}, \ldots, \boldsymbol{\beta}^{(S)}. \]
Example: Simulation exercise
Let’s perform a simulation exercise to evaluate the performance of the previous Algorithm for inference using the Bayesian bootstrap. The data-generating process is defined by two regressors, each distributed as standard normal. The location vector is \(\boldsymbol{\beta} = \left[1 \ 1 \ 1\right]^{\top}\), with a variance of \(\sigma^2 = 1\), and the sample size is 1,000.
The following Algorithm illustrates how to use our GUI to run the Bayesian bootstrap. Our GUI is based on the bayesboot command from the bayesboot package in R. Exercise 11 asks about using this package to perform inference in this simulation and compares the results with those obtained using our GUI with \(S = 10000\).
Algorithm: Bayesian Bootstrap in Linear Regression
Select Univariate Models on the top panel
Select Bootstrap model using the left radio button
Upload the dataset by first selecting whether there is a header in the file, and the kind of separator in the csv file of the dataset (comma, semicolon, or tab). Then, use the Browse button under the Choose File legend. You should see a preview of the dataset
Select number of bootstrap replications using the Range sliders
Select dependent and independent variables using the Formula builder table
Click the Build formula button to generate the formula in R syntax. You can modify the formula in the Main equation box using valid arguments of the formula command structure in R
Click the Go! button
Analyze results
Download posterior chains and diagnostic plots using the Download Posterior Chains and Download Posterior Graphs buttons
The following R code shows how to program the Bayesian bootstrap from scratch. We observe from the results that all 95% credible intervals encompass the population parameters, and the posterior means are close to the population parameters.
rm(list = ls()); set.seed(010101)
#--- Data
N <- 1000; x1 <- runif(N); x2 <- rnorm(N)
X <- cbind(1, x1, x2); B <- c(1, 1, 1); sig2 <- 1
y <- as.numeric(X %*% B + rnorm(N, 0, sqrt(sig2)))
data <- data.frame(y, x1, x2)
#--- Bayesian bootstrap (Rubin) for OLS coefficients
BB <- function(S, df, alpha = 1) {
N <- nrow(df)
Betas <- matrix(NA_real_, nrow = S, ncol = 3)
colnames(Betas) <- c("(Intercept)", "x1", "x2")
pb <- txtProgressBar(min = 1, max = S, style = 3)
for (s in 1:S) {
# One Dirichlet weight vector over the N observations
w <- as.numeric(LaplacesDemon::rdirichlet(1, rep(alpha, N)))
# Weighted least squares = Bayesian bootstrap draw of the OLS functional
fit <- lm(y ~ x1 + x2, data = df, weights = w)
Betas[s, ] <- coef(fit)
setTxtProgressBar(pb, s)
}
close(pb)
return(Betas)
}
S <- 10000
BBs <- BB(S = S, df = data, alpha = 1)## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) 0.9169 0.06359 0.0006359 0.0006359
## x1 1.1738 0.10698 0.0010698 0.0010698
## x2 1.0136 0.03337 0.0003337 0.0003337
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) 0.7934 0.8732 0.917 0.9595 1.042
## x1 0.9647 1.1019 1.173 1.2462 1.381
## x2 0.9477 0.9911 1.014 1.0362 1.079