When trying to detect collinear columns in $X$ a high proportion of cases give a $R_k^2$ close to 1 for independent columns (see figure).
When near multicollinearity arise in a $n\times m$ data matrix $X$, $$X = \left[\begin{array}{ccccc} X_1 & \cdots & X_k & \cdots & X_m \end{array}\right]$$ to be used in the linear regression, $y = X\beta + \varepsilon$, I'm trying to give information to the user on which columns of $X$ are contributing to the trouble.
Say I have identified a near collinear $X$, then my strategy is to normalize $X$ by demeaning the columns and scaling them to have standard deviation $1/\sqrt{n}$. Denoting the normalized matrix by $\tilde{X}$, normalization ensures that $\tilde{X}^\prime \tilde{X}$ is an estimator of the $m\times m$ correlation matrix.
According to section 3.3 in this note from ECB, which is citing Marquardt (1970), the variance inflation factors (VIF's) can be computed by the SVD of $\tilde{X}$, like so,
$$VIF_k = \sum_{j=1}^m \left(\frac{v_{k,j}}{s_j}\right)^2$$
where the SVD of $\tilde{X} = USV^\prime$, with $S=\left(\begin{array}{ccc} s_{1} & & 0\\ & \ddots\\ 0 & & s_{m} \end{array}\right)$ and $V=\left(\begin{array}{ccc} v_{1,1} & \cdots & v_{1,m}\\ \vdots & \ddots & \vdots\\ v_{m,1} & \cdots & v_{m,m} \end{array}\right)$.
The multiple $R^2_k$'s, that is the $R^2$ from regressing $X_k$ on all other columns of $X$ (that is $\{X_{i\neq k}\}$), including an intercept, are given by $$R^2_k=1-\frac{1}{VIF_k}.$$
In Python code:
import numpy as np
def normalize(x):
x_dm = (x - x.mean(axis=0))
return x_dm / np.sqrt((x_dm**2).sum(axis=0))
def vif(right_singular_v, singular_values):
return ((right_singular_v/singular_values)**2).sum(axis=1)
def multiple_r_squared(vif):
return 1 - 1/vif
n = 100
np.random.seed(11)
x1 = np.random.uniform(size=(n, 1))
x2 = np.random.uniform(size=(n, 1))
x3 = np.random.uniform(size=(n, 1))
x4 = np.random.uniform(size=(n, 1))
x5 = 0.5*x3 + 0.5*x4 # x5 is perfectly linear combination of x3 and x4
x = np.concatenate([x1, x2, x3, x4, x5], axis=1)
u, s, v_t = np.linalg.svd(normalize(x), full_matrices=False)
vifs = vif(right_singular_v=v_t.T, singular_values=s)
print(multiple_r_squared(vifs))
>> [ 0.07771458 0.9822463 1. 1. 1. ]
In the above code example x2 is independent of the other x's but the multiple $R^2$ value is calculated at 0.98 which would indicate a that x2 highly correlates with the remaining four variables.
When running a simulation where I repeat the above 10,000 times I get the following distribution of multiple $R^2$'s.
Inspecting the histogram of x2 we see that in approximately 50% of the cases we would reach the false conclusion that x2 is contributing to multicollinearity in x. As the sample size n grows this problem diminishes.
Hoping that my logic and code is correct...
Why is the multiple $R^2$ value of x2 calculated close to 1 in so many cases? Are there, for example, some small sample considerations to take into account, i.e. some asymptotic behavior of the $R_k^2$ estimator that hasn't kicked in yet with only n=100?
