# What exactly is model instability due to multicollinearity?

I've read that model parameters become unstable in case of multicollinearity. Can some one give an example situation of this behavior, and explain why this is happening?

Please use the following multiple linear regression for illustration:

$$y = a_1x_1 + a_2x_2 + b$$

# What Is It?

Here is an example of this behavior. I'm going to write a function to simulate regressions and output their coefficients. We'll look at the coordinate pair of coefficients $$(a_1,a_2)$$ in the case of no collinearity and high collinearity. Here is some code:

    library(tidyverse)

sim <- function(rho){
#Number of samples to draw
N = 50

#Make a covariance matrix
covar = matrix(c(1,rho, rho, 1), byrow = T, nrow = 2)

# Append a column of 1s to N draws from a 2-dimensional
# Gaussian
# With covariance matrix covar
X = cbind(rep(1,N),MASS::mvrnorm(N, mu = c(0,0),
Sigma = covar))

# True betas for our regression
betas = c(1,2,4)

# Make the outcome
y = X%*%betas + rnorm(N,0,1)

# Fit a linear model
model = lm(y ~ X[,2] + X[,3])

# Return a dataframe of the coefficients
return(tibble(a1 = coef(model), a2 = coef(model)))
}

#Run the function 1000 times and stack the results
zero_covar = rerun(1000, sim(0)) %>%
bind_rows

#Same as above, but the covariance in covar matrix
#is now non-zero
high_covar = rerun(1000, sim(0.95)) %>% bind_rows

#plot
zero_covar %>%
ggplot(aes(a1,a2)) +
geom_point(data = high_covar, color = 'red') +
geom_point()


Run that and you get something like This simulation is supposed to simulate the sampling distribution of the coefficients. As we can see, in the case of no collinearity (black dots) the sampling distribution for the coefficients is very tight around the true value of (2,4). The blob is symmetric about this point.

In the case of high collinearity (red dots), the coefficients of the linear model can vary quite a lot! Instability in this case manifests as wildly different coefficient values given the same data generating process.

# Why Is This Happening

Let's take a statistical perspective. The sampling distribution for the coefficients of a linear regression (with enough data) looks like $$\hat{\beta} \sim \mathcal{N}(\beta, \Sigma)$$ The covariance matrix for the above is $$\Sigma = \sigma^{2}\left(X^{\prime} X\right)^{-1}$$ Let's focus for a minute on $$\left(X^{\prime} X\right)$$. If $$X$$ has full rank, then $$\left(X^{\prime} X\right)$$ is a Gram Matrix, which has some special properties. One of those properties is that it has positive eigenvalues. That means we can decompose this matrix product according to eigenvalue decomposition. $$\left(X^{\prime} X\right) = Q\Lambda Q^{-1}$$ Suppose now one of the columns of $$X$$ is highly correlated with another column. Then, one of the eigenvalues of $$X^{\prime}X$$ should be close to 0 (I think). Inverting this product gives us $$\left(X^{\prime} X\right)^{-1} = Q^{-1}\Lambda^{-1} Q$$ Since $$\Lambda$$ is a diagonal matrix, $$\Lambda^{-1}_{jj} = \frac{1}{\Lambda_{jj}}$$. If one of the eigenvalues is really small, then one of the elements of $$\Lambda^{-1}$$ is really big, and so too is the covariance, leading to this instability in the coefficients.

I think I got that right, it has been a long time since I've done linear algebra.

• Looks pretty good, Demetri! But you have activated my inner linear algebra pedant. It's not correct that the Gram matrix has distinct eigenvalues. For example, a full factorial design matrix will have a gram matrix with all identical eigenvalues (and this is not problematic for us). Rather, the full-rank assumption implies strictly positive eigenvalues. (also, the spectral decomp will always exist for a Gram matrix, no need for full-rank or any other assumptions). And you're totally right, collinearity implies vv small eigenvals (close to zero) which cause problems when we invert :) Jul 29 at 17:45
• @JohnMadden Ah see, now I was dubious but made the claim anyway. Let me edit that. Jul 29 at 17:47

Consider an extreme case where the true process is $$y(i)=ax_1(i)+b+e(i)$$, and you plugged a clone variable $$x_2(i)=x_1(i)$$: $$y(i)=a_1x_1(i)+a_2x_2(i)+b+e(i)$$ In this case you see that any pair of coefficients $$a_1,a_2$$ that satisfies $$a_1+a_2=a$$ would be a solution to the problem: $$y(i)=(a_1+a_2)x_1(i)+b+e(i)\Longleftrightarrow y(i)=ax_1(i)+b+e(i)$$

Hence, your regression solver may return you any pair $$(a_1,a-a_1)$$ as a solution. There's your instability. Uncertainty about coefficients is infinite, you can't do any inference about $$a_1$$ or $$a_2$$.

So, in case of perfect multicollinearity in OLS instead of having a single unique solution, you get the continuum of solutions.