What exactly is model instability due to multicollinearity?

Question

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$$

Answer 1

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)[2], a2 = coef(model)[3]))     
}

#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 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.