# How do I compare bootstrapped regression slopes?

Let us assume I have two data sets with n observations of data pairs of independent variable x and dependent variable y each. Let us further assume I want to generate a distribution of regression slopes for each data set by bootstrapping the observations (with replacement) N times and calculating the regression y = a + bx each time. How do I compare the two distributions in order to say the slopes are significantly different? A U-test for testing the difference between the medians of the distributions would be heavily dependent on N, that is, the more often I repeat the bootstrapping the more significant will be the difference. How do I have to calculate the overlap between the distributions to determine a significant difference?

Bootstrapping is done to get a more robust picture of the sampling distribution than that which is assumed by large sample theory. When you bootstrap, there is effectively no limit to the number of bootsamples' you take; in fact you get a better approximation to the sampling distribution the more bootsamples you take. It is common to use $B=10,000$ bootsamples, although there is nothing magical about that number. Furthermore, you don't run a test on the bootsamples; you have an estimate of the sampling distribution--use it directly. Here's an algorithm:

1. take a bootsample of one data set by sampling $n_1$ boot-observations with replacement. [Regarding the comments below, one relevant question is what constitutes a valid 'boot-observation' to use for your bootsample. In fact, there are several legitimate approaches; I will mention two that are robust and allow you to mirror the structure of your data: When you have observational data (i.e., the data were sampled on all dimensions, a boot-observation can be an ordered n-tuple (e.g., a row from your data set). For example, if you have one predictor variable and one response variable, you would sample $n_1$ $(x,y)$ ordered pairs. On the other hand, when working with experimental data, predictor variable values were not sampled, but experimental units were assigned to intended levels of each predictor variable. In a case like this, you can sample $n_{1j}$ $y$ values from within each of the $j$ levels of your predictor variable, then pair those $y$s with the corresponding value of that predictor level. In this manner, you would not sample over $X$.]
2. fit your regression model and store the slope estimate (call it $\hat\beta_1$)
3. take a bootsample of the other data set by sampling $n_2$ boot-observations with replacement
4. fit the other regression model and store the slope estimate (call it $\hat\beta_2$)
5. form a statistic from the two estimates (suggestion: use the slope difference $\hat\beta_1-\hat\beta_2$)
6. store the statistic and dump the other info so as not to waste memory
7. repeat steps 1 - 6, $B=10,000$ times
8. sort the bootstrapped sampling distribution of slope differences
9. compute the % of the bsd that overlaps 0 (whichever is smaller, the right tail % or the left tail %)
10. multiply this percentage by 2

The logic of this algorithm as a statistical test is fundamentally similar to classical tests (e.g., t-tests) but you are not assuming the the data or the resulting sampling distributions have any particular distribution. (For example, you are not assuming normality.) The primary assumption you are making is that your data are representative of the population you sampled from / want to generalize to. That is, the sample distribution is similar to the population distribution. Note that, if your data are not related to the population you're interested in, you are flat out of luck.

Some people worry about using, e.g., a regression model to determine the slope if you're not willing to assume normality. However, this concern is mistaken. The Gauss-Markov theorem tells us that the estimate is unbiased (i.e., centered on the true value), so it's fine. The lack of normality simply means that the true sampling distribution may be different from the theoretically posited one, and so the p-values are invalid. The bootstrapping procedure gives you a way to deal with this issue.

Two other issues regarding bootstrapping: If the classical assumptions are met, bootstrapping is less efficient (i.e., has less power) than a parametric test. Second, bootstrapping works best when you are exploring near the center of a distribution: means and medians are good, quartiles not so good, bootstrapping the min or max necessarily fail. Regarding the first point, you may not need to bootstrap in your situation; regarding the second point, bootstrapping the slope is perfectly fine.

• Although I may well be wrong, I thought the bootstrap in regression had to be on the residuals rather than on the raw data, to be validated... Nov 22 '11 at 14:05
• @Xi'an, I've been wrong before myself, but I don't understand why you think only bootstrapping residuals is valid. Efron & Tibshirani (1994) section 9.5 says "Bootstrapping pairs is less sensitive to assumptions than bootstrapping residuals. The standard error obtained by bootstrapping pairs gives reasonable answers even if [the probability structure of the linear model] is completely wrong." The implication is that application of the bootstrap is more robust, although they imply it can be less efficient in some cases. Nov 25 '11 at 2:46
• My worry with bootstrapping the pairs is that you also incorporate the distribution of the predictors, which is usually left outside of the picture in regular linear models. That's why I always teach my students to bootstrap only the residuals. Nov 25 '11 at 8:11
• @Xi'an, that is a reasonable point, I suppose I was assuming an observational data structure. I have edited my answer to add more detail about these concerns. However, I do not see how that implies that bootstrapping pairs is necessarily invalid. Nov 25 '11 at 17:23
• Forcing a pairing between two independent sets of data is artificial and inefficient. You can do much better than that!
– whuber
Nov 25 '11 at 17:57

You can combine the two data sets into one regression. Let $s_i$ be an indicator for being in the first data set. Then run the regression $$\begin{equation*}y_i = \beta_0 + \beta_1 x_i + \beta_2 s_i + \beta_3 s_i x_i + \epsilon_i \end{equation*}$$ Note that the interpretation of $\beta_3$ is the difference in slopes from the separate regressions: \begin{align*} \text{E}[y_i \mid x, s_i = 1] &= (\beta_0 + \beta_2) + (\beta_1 + \beta_3) x_i \\ \text{E}[y_i \mid x, s_i = 0] &= \beta_0 + \beta_1 x_i. \end{align*} You can bootstrap the distribution of $\beta_3$ if you want or just use standard testing procedures (normal/t). If using analytical solutions, you need to either assume homoskedasticity across groups or correct for heteroskedasticity. For bootstrapping to be robust to this, you need to choose $n$ observations randomly among the first group and $n$ among the second, rather than $2n$ from the whole population.

If you have correlation among the error terms, you may need to alter this procedure a bit, so write back if that is the case.

You can generalize this approach to the seemingly unrelated regressions (SUR) framework. This approach still allows the coefficients for the intercept and the slope to be arbitrarily different in the two data sets.

• It's a good idea. But doesn't this also assume the two regressions have iid errors?
– whuber
Nov 21 '11 at 15:48
• Good point. It requires that there not be different variances for the errors by group and that the errors not be correlated in the different groups. Nov 21 '11 at 17:23

Doing everything in one regression is neat, and the assumption of independence is important. But calculating the point estimates in this way does not require constant variance. Try this R code;

x <- rbinom(100, 1, 0.5)
z <- rnorm(100)
y <- rnorm(100)
coef(lm(y~x*z))
coef(lm(y~z, subset= x==1))[1] - coef(lm(y~z, subset= x==0))[1]
coef(lm(y~z, subset= x==1))[2] - coef(lm(y~z, subset= x==0))[2]
`

We get the same point estimate either way. Estimates of standard error may require constant variance (depending on which one you use) but the bootstrapping considered here doesn't use estimated standard errors.

• If you're going to test whether the difference of slopes is zero (as in @Charlie's reply, to which you seem to be following up), you need an accurate, valid estimate of standard errors. It doesn't matter whether you bootstrap that estimate or otherwise.
– whuber
Nov 22 '11 at 3:45