Why doesn't collinearity affect the predictions?

I have read in many places that collinearity doesn't affect the predictions. It only affects the coefficient tests and confidence interval. As a result it cannot be used for causal inference but for making predictions it is 'safe' to use. Here there is a proof for the case when the predictors have perfect correlation. But why is this the case for imperfectly correlated predictors?

Let's say we have two predictors with correlation $$0.8$$, how does inclusion of these two in a multiple regression affect the predictions? Below I have an example with $$0.91$$ correlation between two predictors, and it looks like there is some difference between predictions, especially larger on test data. But is it possible to have a formal analysis (like the one for perfect correlation in the link above).

> x1 <- 1:100 + rnorm(100)
> x2 <- x1^3 - rnorm(100, 50, 50)
> cor(x1, x2)
[1] 0.9179052
> y <- 201:300 + rnorm(100, 0, 20)
> cor(y, x1)
[1] 0.8463703
> cor(y, x2)
[1] 0.7621107
> d <- data.frame(x1= x1, x2= x2, y= y)
> m <- lm(y~x1, data = d[1:50,])
> m2 <- lm(y~x2, data = d[1:50,])
> m3 <- lm(y~x1+x2, data = d[1:50,])

> #on training data
> p1 <- predict(m, d[1,50])
> p2 <- predict(m2, d[1,50])
> p3 <- predict(m3, d[1,50])
> summary(p1 - p2)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-10.759  -3.269   1.320   0.000   3.872   5.166
> summary(p1 - p3)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-0.8312 -0.6590 -0.2200  0.0000  0.5984  1.5835
> summary(p2 - p3)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-5.979  -4.550  -1.755   0.000   3.877  12.052

> #on test data
> p4 <- predict(m, d[51:100,])
> p5 <- predict(m2, d[51:100,])
> p6 <- predict(m3, d[51:100,])
> summary(p4 - p5)
Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
-232.308 -137.667  -76.297  -90.126  -33.208   -7.692
> summary(p5 - p6)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
9.709  40.101  91.202 107.497 163.833 275.705
> summary(p4 - p6)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
2.016   6.893  14.905  17.371  26.166  43.397

> #model summaries
> summary(m)

Call:
lm(formula = y ~ x1, data = d[1:50, ])

Residuals:
Min      1Q  Median      3Q     Max
-45.424 -10.758   0.179  13.969  31.637

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 204.0189     5.5249  36.927  < 2e-16 ***
x1            0.8227     0.1898   4.335 7.43e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19.04 on 48 degrees of freedom
Multiple R-squared:  0.2813,    Adjusted R-squared:  0.2664
F-statistic: 18.79 on 1 and 48 DF,  p-value: 7.435e-05

> summary(m2)

Call:
lm(formula = y ~ x2, data = d[1:50, ])

Residuals:
Min      1Q  Median      3Q     Max
-44.193 -14.038   0.306  15.814  35.632

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.157e+02  3.736e+00  57.725  < 2e-16 ***
x2          2.923e-04  7.808e-05   3.743 0.000486 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19.76 on 48 degrees of freedom
Multiple R-squared:  0.226, Adjusted R-squared:  0.2098
F-statistic: 14.01 on 1 and 48 DF,  p-value: 0.0004857

> summary(m3)

Call:
lm(formula = y ~ x1 + x2, data = d[1:50, ])

Residuals:
Min      1Q  Median      3Q     Max
-45.848 -11.305   0.024  13.404  32.701

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.026e+02  7.700e+00  26.310   <2e-16 ***
x1           9.439e-01  4.908e-01   1.923   0.0605 .
x2          -5.214e-05  1.945e-04  -0.268   0.7898
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19.23 on 47 degrees of freedom
Multiple R-squared:  0.2824,    Adjusted R-squared:  0.2519
F-statistic: 9.249 on 2 and 47 DF,  p-value: 0.0004101


Second question is how would inclusion of new variables that are linear combination of these two correlated variable affect the predictions? Let's say doing a regression on $$x1$$,$$x2$$, and $$x3$$ where $$x3=x2+x1$$.

> d <- cbind(d, data.frame(x3 = d$$x1 + d$$x2, x4 = d$$x1 - d$$x2))
> m4 <- lm(y~x3+x4, data = d[1:50,])
> p7 <- predict(m4, d[1:50,])
> p8 <- predict(m4, d[51:100,])
> summary(p7-p3)
Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-1.623e-11 -9.898e-12 -4.107e-12  3.604e-13  8.413e-12  3.175e-11
> summary(p8-p6)
Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-1.003e-11  3.312e-11  8.608e-11  1.299e-10  2.283e-10  3.933e-10

• Your second paragraph seems unlikely to me unless something special is going on. Unless there is perfect collinearity of some sort (including more predictors than observations), I would have thought changing the coefficients should change some of the predictions – Henry Feb 3 '20 at 17:43
• I think you have in mind that the term "collinear" refers to correlated predictors even if the correlation is not 1. Chances are that what you read about "doesn't affect the predictions" refers to perfect collinearity, i.e., having "perfect" linear dependency between predictors. – Lewian Feb 3 '20 at 18:14
• Henry, seems like you were right, maybe there was something special in my example before, I re-did the analysis and now there is a difference in predictions, more so on test data rather than training data. (I guess I can see why, because my test data are "extrapolations" having x2 ~ x1^3 and test data on the second half of the data). I have changed the question a bit now. – armen Feb 3 '20 at 18:17
• Also "collinearity doesn't affect the predictions" is somewhat ambiguous. It may just mean that under (perfect) collinearity you can still properly compute unique predictions, which you surely can do if you don't have perfect collinearity (that case would be a triviality if "affect predictions" is interpreted in this way). You seem to interpret the statement in terms of whether inclusion or exclusion of variables affects prediction under collinearity. Inclusion of new variables that are perfect linear combinations of existing ones won't change predictions. Otherwise it may well. – Lewian Feb 3 '20 at 18:21
• @armen: Because if you are adding a variable that is a linear combination of others, you are not adding new information. The information in the x from which to predict the y remains the same. – Lewian Feb 3 '20 at 18:40

It depends on how those predictions (whatfor) are made.

Say you do a study on the influence of alcoholic products on health. Say you measure the number of occurances of cardiovascular disease as dependent variable and you have, among other beverages, separately white wine consumption and red wine consumption as independent variables.

For a lot of people, the case is that they drink both if they drink any (the consumption of red and white wine correlates). If there is an effect of wine then it is a bit interchangeable to put this effect on white or red. So there might be in some studies (if they have red and white wine correlated) not much clarity about whether it is red or white.

However, on the one hand, for most predictions it doesn't really matter. If you would erroneously put the effect on the one variable then you are still predicting indirectly ok because this wrong variable correlates with the correct variable.

On the other hand, if you are making the predictions for cases where this correlation is absent (e.g. a group/person that is drinking only white wine or only red wine) then you may predict badly. The estimate/prediction may have bad external validity.

See also the example in my answer to Does LASSO suffer from the same problems stepwise regression does?

It shows fitting of the curve allong with the estimated coefficients. Some of those coefficients are estimated badly but because the wrongly estimated components are so close (correlating) to the correct components this matters very little for the fit of the curve.

How does inclusion of these two in a multiple regression affect the predictions

In this question Is ridge regression useless in high dimensions ($n \ll p$)? How can OLS fail to overfit? you can read how additional parameters work as some sort of additional noise, but also can actually work regularizing.

This is a situation where mathematical abstraction is a huge help. In the following I will not introduce any new ideas, nor carry out any calculations, but will only exploit basic, simple definitions of linear algebra to present an effective way of thinking about collinearity.

The definitions needed to understand this are vector space, subspace, linear combination, linear transformation, kernel, and orthogonal projection.

One abstraction that works well is to view the the columns of the $$n\times p$$ design matrix $$X$$ as elements of the vector space $$\mathbb{R}^n.$$ For any $$p$$-vector $$\beta,$$ interpret the expression

$$X\beta\tag{*}$$

as a linear combination of these columns. The set of all such possible linear combinations is a subspace $$\mathbb{R}[X] \subset \mathbb{R}^n.$$ (Its dimension is at most $$p$$ but could be smaller. We say $$X$$ is "collinear" when $$\operatorname{dim}(\mathbb{R}[X])$$ is strictly less than $$p.$$)

The expression $$(*)$$ defines a linear transformation, which to avoid abusing notation I will give a new name,

$$\mathcal{L}_X: \mathbb{R}^p \to \mathbb{R}[X] \subset \mathbb{R}^n;\quad \mathcal{L}_X(\beta) = X\beta.$$

Let's forget about this transformation a moment in order to focus on the image of $$\mathcal{L}_X$$ in $$\mathbb{R}^n.$$ To do so, call this subspace $$V = \mathbb{R}[X]$$ to emphasize that it's just some definite subspace and to de-emphasize that $$X$$ was used to figure exactly which subspace $$V$$ is. The response vector $$y$$ can also be viewed as an element of $$\mathbb{R}^n.$$

I hope the following characterization requires no proof (because it merely restates well-known results):

Ordinary Least Squares regression seeks a $$\hat y \in V$$ of minimal distance to $$y.$$ The solution is unique and is found by projecting $$y$$ orthogonally onto $$V.$$

Let us call this projection $$\operatorname{proj}:\mathbb{R}^n \to V,$$ so that

$$\hat y = \operatorname{proj}_V(y).$$

Note that $$\operatorname{proj}$$ depends only on (1) the metric on $$\mathbb{R}^n,$$ up to a scalar multiple (this is given by some assumption about the covariances of the error terms in the model) and (2) the subspace $$V.$$

This figure sketches the vector spaces and transformations among them. The model fitting occurs in $$\mathbb{R}^n$$ shown at the upper right, whereas how it is described depends on the objects at the lower left. The vectors $$x_i\in\mathbb{R}^p$$ are mapped to the corresponding columns of $$X,$$ which are elements of $$V\subset \mathbb{R}^n.$$ $$\hat\beta$$ is mapped via $$\mathcal{L}_X$$ to $$\hat y.$$

The point of these definitional trivialities is to make it clear that this characterization depends on $$X$$ only through the subspace $$V=\mathbb{R}[X]$$ it determines. Some immediate consequences are:

1. The solution does not depend on whether the columns of $$X$$ are orthogonal.

2. Introducing additional columns into $$X$$ that are already elements of $$V$$ changes nothing.

3. Parameter estimates $$\hat\beta$$ are obtained by pulling $$\hat y$$ back via $$\mathcal{L}_X.$$ Equivalently, all such $$\hat\beta$$ satisfy the equation $$\mathcal{L}_X(\hat\beta) = \operatorname{proj}_V(y) = \hat y.$$

4. Any two such parameter estimates $$\hat\beta_1, \hat\beta_2\in\mathbb{R}^p$$ must therefore differ by an element of the kernel of $$\operatorname{proj}_V,$$ because $$\operatorname{proj}_V(\hat\beta_2 - \hat\beta_1) = \operatorname{proj}_V(\hat\beta_2) - \operatorname{proj}_V(\hat\beta_1) = \hat y - \hat y = 0 \in \mathbb{R}^n.$$

Point (3) shows the prediction $$\hat y$$ does not depend on $$\hat\beta$$ and point (4) characterizes the extent to which $$\hat\beta$$ may be incompletely characterized. Of course, when $$\operatorname{dim} V = p,$$ $$\operatorname{ker}(\operatorname{proj}_V) = 0$$ and therefore $$\hat\beta$$ is uniquely determined.

• "this is given by some assumption about the covariances of the error terms in the model" would you mind expanding on this connection? you mean if the covariance structure of the error terms were something other than uncorrelated, you would use some mahalonobis-type distance rather than the standard metric? – Hasse1987 Feb 4 '20 at 22:08
• @Hasse Yes--but nothing changes conceptually. All metrics on $\mathbb{R}^n$ given by non-degenerate, positive-definite quadratic forms are equivalent: you may view them as differing only in how one describes the usual Euclidean metric. (See stats.stackexchange.com/a/62147/919 for an explanation of this equivalence, which is carried out in a setting where the covariance is determined by a dataset--but see the last technical comment at the very end.) The covariances of the error terms determine this metric. – whuber Feb 4 '20 at 22:20
• Thank you for this explanation, i really appreciated the approach but was a bit lost at the end ( not so trivial for me it seems). What do you means by the vectors xi are mapped to the corresponding columns of X (what are these xi by the way). Is it also possible to elaborate on (2): I mean if we add a column in X which is partially colinear with other columns it will definetly change V and therefore yhat (3): How does this point show that prediction does not depends on Beta estimate (very contrintuitive from the regression equation, if i change beta hat it will change yhat) . Best Alexis – Alexis Mar 22 at 1:46
• @Alexis I appreciate your interest, but anticipated this answer would not be suitable for everyone. As announced at the beginning, once you learn about the concepts of "vector space, subspace, linear combination, linear transformation, kernel, and orthogonal projection," this should become accessible to you. – whuber Mar 22 at 15:11