# Picking random regression lines from the distribution of OLS regressions

I regressed a series of data with OLS, and I calculated the standard error (SE) of the slope and of the intercept. Moreover, I calculated the confidence bands for the regression following the equations in the answer to this question. Now I am wondering how I can generate a set of random regression lines within this distribution of regression lines that is illustrated by the 95% confidence interval.

In order to illustrate what I mean, here is a short summary of my problem. I have a series of data (dataset A; blue points in below figure) that deviate from a linear trend defined by another dataset (dataset B; red points in below figure). Using some constraints from chemistry, I can correct A by projecting each point onto the linear trend defined by B (red regression line in below figure). In order to consider the uncertainty in the regression of dataset B in the corrected values for dataset A, I was thinking of picking lines at random from the distribution of regressions of dataset B, and projecting A onto each of these lines. I just do not know how to pick lines that are representative of the distribution of regressions.

Below is now a zoom into the same plot as above. Now, the black lines show 100 regression lines picked at random. These lines were generated by choosing random values for the slope and intercept from normal distributions with mean and standard errors obtained from the regression.

Here is the same with 1000 lines

There seem to be more than 5% of the lines that fall out of the 95% confidence bands, and the method of picking a random slope and a random intercept value does not work. This makes sense, considering that slope and intercept are not independent. Is there a simple way to generate regressions that are representative of the true distribution of regression lines around the "best fit"?

Thanks and please let me know if any of this is unclear

• What is a 'mixing line'? Can you say more about these datasets? – mkt - Reinstate Monica Jul 22 '19 at 9:06
• O.K. - I see how that is confusing. For the purpose of this question, it is just a regression line that I project dataset A onto. I will edit that in the text – NAB Jul 22 '19 at 9:17

In a linear regression model, the joint sampling distribution of the model parameters is described in this post and this one. Specifically, $$\hat{\beta}\;|\;X\sim \mathrm{N}\left(\beta, \sigma^{2}(X^TX)^{-1}\right)$$ where $$X$$ is the design matrix. In a simple linear regression with only one predictor and an intercept, the design matrix has dimensions $$N\times 2$$ (i.e. it has two columns and $$N$$ rows) where $$N$$ is the sample size. You could then sample from a bivariate normal distribution to generate regression lines. Another possibility would be to bootstrap.

If the predictor variable was centered on its sample mean before the regression, then the covariance between the intercept and the slope is $$0$$ (for more information, look here). More general: If the sample mean of the predictor variable $$x$$ is not $$0$$ ($$\bar{x}\neq 0$$), the intercept and the slope are correlated. If the predictor variable was not centered, the intercept and the slope are not independent. The following graphic depicts the 80%, 95% and 99% confidence ellipse for the slope and the intercept in the case that the predictor variable was not centered:

If the predictor were centered, the corresponding confidence ellipses look different:

Here are two graphics, the first was generated by sampling 500 times from the bivariate normal distribution and the second one was generated by 500 nonparametric bootstrap replicates. The red line denotes the ols fit and the dashed black lines denote the 95% confidence intervals. The individual resampled regression lines are drawn in a translucent black so that the color is darker in regions were they overlap. You can see clearly that the majority of the resampled regression lines lie within the 95% confidence band, as we would expect.

Here is the R code I used:

library(MASS)
library(scales)
library(ggplot2)
library(reshape2)

mod <- lm(Infant.Mortality~Fertility, data = swiss)

Sigma <- vcov(mod) # The covariance matrix of the estimates
mu <- coef(mod) # The estimates

set.seed(142857)

n_sim <- 500      # Number of replicates

x <- mvrnorm(n = n_sim, mu = mu, Sigma = Sigma) # Sample from the bivariate normal distribution

rbind(x, c(mu))

plot_frame <- data.frame(
int = x[, 1]
, slope = x[, 2]
)

pred_frame <- data.frame(
Fertility = seq(min(swiss$$Fertility), max(swiss$$Fertility), length.out = 1000)
)

pred_frame$$lwr <- predict(mod, newdata = pred_frame, interval = "confidence")[, 2] pred_frame$$upr <- predict(mod, newdata = pred_frame, interval = "confidence")[, 3]

theme_set(theme_bw())
p <- ggplot(data = swiss, aes(x = Fertility, y = Infant.Mortality)) +
geom_point(size = 3, colour = "#5CC3F0", alpha = 0.5, pch = 16) +
geom_abline(data = plot_frame, aes(intercept = int, slope = slope), col = "black", alpha = 0.05) +
geom_line(data = pred_frame, aes(x = Fertility, y = lwr), linetype = 2, size = 1.7) +
geom_line(data = pred_frame, aes(x = Fertility, y = upr), linetype = 2, size = 1.7) +
stat_smooth(method = "lm", se = FALSE, col = "red", size = 1.5, fullrange = TRUE) +
xlab("Standardized fertility measure") +
ylab("Live births who live less than 1 year") +
ggtitle("Bivariate normal distribution") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
theme(
axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
axis.title.x=element_text(colour = "black", size = 17),
axis.text.x=element_text(colour = "black", size=15),
axis.text.y=element_text(colour = "black", size=15),
legend.position="right",
legend.text=element_text(size=12.5),
legend.key=element_blank(),
plot.title = element_text(face = "bold"),
strip.text.x=element_text(size=15)
)

#===============================================================
# Bootstrap
#===============================================================

set.seed(142857)

n_boot <- 500 # Number of bootstrap replicates

boot_res <- matrix(NA, ncol = 2, nrow = n_boot) # Matrix to store the bootstrap coefficients

# Start the bootstrap

for (i in seq_len(n_boot)) {
boot_res[i, ] <- coef(lm(Infant.Mortality~Fertility, data = swiss[sample(seq_len(dim(swiss)[1]), replace = TRUE), ]))
}

boot_res <- data.frame(boot_res)
names(boot_res) <- c("int", "slope")

theme_set(theme_bw())
p <- ggplot(data = swiss, aes(x = Fertility, y = Infant.Mortality)) +
geom_point(size = 3, colour = "#5CC3F0", alpha = 0.5, pch = 16) +
geom_abline(data = boot_res, aes(intercept = int, slope = slope), col = "black", alpha = 0.05) +
geom_line(data = pred_frame, aes(x = Fertility, y = lwr), linetype = 2, size = 1.7) +
geom_line(data = pred_frame, aes(x = Fertility, y = upr), linetype = 2, size = 1.7) +
stat_smooth(method = "lm", se = FALSE, col = "red", size = 1.5, fullrange = TRUE) +
xlab("Standardized fertility measure") +
ylab("Live births who live less than 1 year") +
ggtitle("Bootstrap") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
theme(
axis.title.y=element_text(colour = "black", size = 17, hjust = 0.5, margin=margin(0,12,0,0)),
axis.title.x=element_text(colour = "black", size = 17),
axis.text.x=element_text(colour = "black", size=15),
axis.text.y=element_text(colour = "black", size=15),
legend.position="right",
legend.text=element_text(size=12.5),
legend.key=element_blank(),
plot.title = element_text(face = "bold"),
strip.text.x=element_text(size=15)
)

• This is fantastic and looks exactly like I was looking for! It will take me a moment to go through the answer and code, and I will see tomorrow if I can make it work! Thank you so much! – NAB Jul 22 '19 at 19:19
• I had a look at it and it worked beautifully! Just one weird confusion: I use matlab and, there is a function in matlab: [beta,Sigma,E,CovB,logL] = mvregress(X,Y), where CovB is supposedly the "variance-covariance matrix of the regression coefficients" according to the documentation: de.mathworks.com/help/stats/mvregress.html With that function, I get the same coefficients (beta), but not the same covariance matrix than vcov (and vcov gives the same matrix than one I calculated manually using matrix algebra). Are there different definitions of the "variance-covariance matrix"? – NAB Jul 23 '19 at 8:08
• @NAB No, the variance-covariance matrix should be the same. For a simple linear regression, the formula for the covariance matrix is actually quite simple, look here. You could cross-check by calculating the covariance matrix "by hand", so to speak. I don't have access to Matlab, but the help page you linked to is for multivariate regression (i.e. multiple dependent variables). I'm not sure if that could be the source of the difference. But I have found this tutorial that might be helpful. – COOLSerdash Jul 23 '19 at 8:14
• Yeah - I calculated it by hand and got the same result as the vcor function in R (I never really used R, but your code was well-explained) - s0, just odd to have that Matlab function give a different result. But I suppose that is a different topic; I might just be missing something. In any case, thanks for the help! – NAB Jul 23 '19 at 10:01
• @NAB No problem. – COOLSerdash Jul 23 '19 at 10:24