Practical limits to collinearity problems?

Question

Collinear independent variables can have undesirable effects on the interpretation of coefficients in a linear model. Indeed, for two perfectly-correlated predictors, the coefficients are not uniquely determined, leaving a single degree of freedom by which they may vary.

However, adding in even a little bit of randomness gives solutions to the OLS equations that are not too far from assigning the correlated independent variables equal weight.

An example. The following (python) code creates a simple linear relationship between a single independent variable and a single dependent variable. It then copies the independent variable to generate a second independent variable and adds a tiny amount of independent gaussian noise to both. Thus, the second variable is nearly-perfectly correlated with the first variable. When we do this a number of times and plot the resulting coefficients, they tend to center around 2.5 and 2.5:

import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model

coefs = []                                # to hold the coefficient of all of the OLS fits
for i in range(1000):                     # run OLS a bunch to see what the coefficients do
    X = np.linspace(0,1,100)              # some independent variable
    y = 5*X + 4                           # some dependent variable
    X2 = X                                # a new independent variable that is correlated with the first
    XX = (                                
        np.stack((X,X2)).T +              # stack the independent variables
        np.random.normal(0,0.01,(100,2)   # add noise
    )
    lr = sklearn.linear_model.LinearRegression()
    lr.fit(XX,y)
    coefs.append(lr.coef_)                # get the coefficients of an OLS linear regression

coefs = np.array(coefs)
plt.scatter(coefs[:,0], coefs[:,1])
plt.show()

Why am I relatively unlikely to get, say, an X1 coefficient of -105 and an X2 coefficient of 110? Those add up to 5, but there is something pushing the results toward 2.5, 2.5. What is the intuition behind this phenomenon, and what practical implications does this have when faced with collinear independent variables?

UPDATE:

Altering the amount of Gaussian noise to have a much smaller variance ($\sigma = 0.00000001$ as opposed to $\sigma=0.01$ in the example above):

And $\sigma=0.0000000000000001$:

UPDATEUPDATE:

Interest in how the amount of added noise plays into this phenomenon encouraged me to make the following graph. I ran the above experiment for various quantities of added noise, from $2^{1}$ to $2^{-60}$ on an exponential scale. I then measured the std of the X1 coefficient in the above experiment for each of those trials. Results are given here (should read $-\log_2$ on the x-axis):

I suspect that at very very small values of added noise, there are numerical problems due to limits of the computer architecture. Accounting for that, it seems from this graph that the distribution of OLS coefficients for X1 and X2 converges to a Gaussian distribution with mean 2.5 and std 0.25 or so. Thus, this seems like "Why not coefficients (-105,110)?" is still open.

Zooming in on the weird behavior for very small values:

Which seems like it's doing a sane thing. That doesn't mean that it's not floating-point errors, but it's not clear that that is what's happening.

Answer 1

Why am I relatively unlikely to get, say, an x_1 coefficient of -105 and an x_2 coefficient of 110? Those add up to 5, but there is something pushing the results toward 2.5, 2.5.

Linear combinations of your two variables $x_1 = x + \epsilon_1$ and $x_2 = x+\epsilon_2$ can be described like:

$$ \frac{a+b}{2} (x + \epsilon_1) + \frac{a-b}{2}(x + \epsilon_2) = a x + \frac{1}{2} a (\epsilon_1 + \epsilon_2) + \frac{1}{2} b (\epsilon_1 - \epsilon_2)$$

• The parameter a will be approximately equal to the parameter associated with the variale $x$. In your case that is $a = 5$.

• The parameter b will be related to the variance and correlation of $y=(\epsilon_1 + \epsilon_2)$ and $z=(\epsilon_1 - \epsilon_2)$ by: $$\text{Var} \left( a y + b z \right) = a^2 \text{Var} (y) + b^2 \text{Var}(z) + 2ab \sqrt{\text{Var}(y)\text{Var}(z)} \rho_{y,z} $$ note that $y=(\epsilon_1 + \epsilon_2)$ and $z=(\epsilon_1 - \epsilon_2)$ are iid distribtued variables and $\rho_{y,z}$ will be distributed around zero.

  So mostly you will be close to $b=0$

So in the case of -105 and 110 you would get larger contributions from the error terms, which only gets 'undone' when there is an a strong correlation in the particular sample of the error terms.

Influence of $\sigma$

I can model as well the influence of $\sigma$ but I do not get the same pattern as you. Below you see that with larger variance, the sum of the parameters will be smaller than $5$ (to decrease the effect of the error terms) and also the difference will be smaller as it relates to the size of the sum of the parameters. But, I do not see why the parameters would go from 0 to 5 as in your last graph.

---

This is for thousand repetitions of the data:

$$\begin{array}{rcl} x_1 &=& x + \epsilon_1 \\ x_2 &=& x + \epsilon_2 \\ y &=& 5x + 4 \end{array}$$

where $x$ is a vector of size $n=100$ varying from 0 to 1, $\epsilon_1$ and $\epsilon_2$ are Gaussian random noise.

Which is modeled as a linear model minimizing the sum of least squared error $\epsilon$

$$y = a x_1 + b x_2 + \epsilon$$

Answer 2

In the way you've added noise, you could write $x' = x +\epsilon$ (where $\epsilon$ is a normally distributed variable representing noise). Furthermore, $y=5x+4+\eta$, where $\eta$ is another Gaussian noise term.

You are trying to fit a regression of the form $Ax + Bx' +C$, and you know that you're fitting it to a target variable which was generated by $y=5x+4+\eta$, thus :

$Ax + B(x + \epsilon) + C= 5x+4+\eta$

or

$(A+B)x + C + B\epsilon = 5x + 4 + \eta$

Hopefully, this sheds some light on what's happening (together with your numerics), even if it doesn't formally show it. If $\epsilon$ is so small that $B\epsilon$ is also very small, then pretty much any combination of (A,B) which satisfies $A+B=5$ will likely be a good fit, and which one is best is a question of noise. The large $\epsilon$ is, the more $B\epsilon$ will in general make the LHS fluctuate wrt to the RHS and the more B will go to zero ( a trivial restatement of the fact that x is the truly explanatory variable and B is only correlated, and as you shrink that correlation, the regression will find it easier to learn this)

It's less clear to me what happens when $\epsilon$ and $\eta$ are of similar sizes. In general there will be competing effects where you might be able to match the data better with non-zero B. What happens if you simulate a larger data set? Does the effect still hold?