Variance inflation factors not equal to $C_{jj}$

Question

I am following a book which states that the diagonal elements of $C = (X'X)^{-1}$ are called the variance inflation factors:

$$ VIF_j = C_{jj} = \frac{1}{1-R^2_j} $$

where $R^2_j$ is the coefficient of determination if the $x_j$ variable is regressed on the rest of the variables. When I calculate this $C$ matrix, I get different numbers in the diagonals versus the numbers returned by car::vif(m) which also match the VIFs in the book.

I can give you a dummy example, but here is a concrete example from the book:

y = c(49.0,50.2,50.5,48.5,47.5,44.5,28.0,31.5,34.5,35.0,38.0,
                       38.5,15.0,17.0,20.5,29.5)
x1 = c(1300,1300,1300,1300,1300,1300,1200,1200,1200,1200,1200,1200,
                 1100,1100,1100,1100)
x2 = c(7.5,9.0,11.0,13.5,17.0,23.0,5.3,7.5,11.0,13.5,17.0,23.0,5.3,7.5,11.0,17.0)
x3 = c(0.0120,0.0120,0.0115,0.0130,0.0135,0.0120,0.0400,0.0380,0.0320,
                 0.0260,0.0340,0.0410,0.0840,0.0980,0.0920,0.0860)

# standardize variables - zero mean, unit variance
x1 = scale(x1)
x2 = scale(x2)
x3 = scale(x3)

x1x2 = x1*x2
x1x3 = x1*x3
x2x3 = x2*x3

x12 = x1^2
x22 = x2^2
x32 = x3^2

m = lm(y ~ x1 + x2 + x3 + x1x2 + x1x3 + x2x3 + x12 + x22 + x32)

X = cbind(x1, x2, x3, x12, x22, x32, x1x2, x1x3, x2x3)
C = solve(t(X)%*%X)
diag(C)

The last line gives me:

          x1           x2           x3          x12          x22          x32         x1x2 
 18.23702168   0.10311100  33.05102775 181.63651371   0.09737908  72.77598020   2.62389599 
        x1x3         x2x3 
527.66678556   3.28427247 

However, calling car::vif(m) gives me the following which the book also states as the correct VIFs:

         x1          x2          x3        x1x2        x1x3        x2x3         x12         x22 
 375.247759    1.740631  680.280039   31.037059 6563.345193   35.611286 1762.575365    3.164318 
        x32 
1156.766284

I don't quite understand what I'm doing wrong.

I can see that if I regress $x_1$ against all other regressors, I get $R^2 = 0.9973351$ and then $1/(1-R^2) = 375.2486$ which is indeed the VIF from the vif() function. But this is different to the $C_{jj}$ coefficient.

Answer 1

Your matrix X is not the correct design matrix of the model, because it should contain the column of one's corresponding to the intercept. To get the design matrix, you can do X <- model.matrix(m).

Answer 2

As stated in this CV post:

The standardized interaction term should be the standardized version of the product of the two original variables, not the product of the two standardized variables.

Also, the scaling should be $x_{ij}\mapsto (x_{ij}-\bar x_j)/\sqrt{\sum_i (x_{ij}-\bar x_j)^2}$ (cf. $\rm [I]$).

x1s <- scale(x1, scale = sqrt(var(x1)*(length(x1)-1)))
x2s <- scale(x2, scale = sqrt(var(x2)*(length(x2)-1)))
x3s <- scale(x3, scale = sqrt(var(x3)*(length(x3)-1)))

x1x2s <- scale(x1*x2, scale = sqrt(var(x1*x2)*(length(x1)-1)))
x1x3s <- scale(x1*x3, scale = sqrt(var(x1*x3)*(length(x1)-1)))
x2x3s <- scale(x2*x3, scale = sqrt(var(x2*x3)*(length(x1)-1)))

x12s <- scale(x1^2, scale = sqrt(var(x1^2)*(length(x1)-1)))
x22s <- scale(x2^2, scale = sqrt(var(x2^2)*(length(x2)-1)))
x32s <- scale(x3^2, scale = sqrt(var(x3^2)*(length(x3)-1)))

ms <- lm(y ~ x1s + x2s + x3s + x1x2s + x1x3s + x2x3s + x12s + x22s + x32s)

Xs <- cbind(x1s, x2s, x3s, x1x2s, x1x3s, x2x3s, x12s, x22s, x32s)

Cs <- solve(t(Xs)%*%Xs)

diag(Cs)

[1] 2.856749e+06 1.095614e+04 2.017163e+06 9.802903e+03 1.428092e+06 2.403594e+02
[7] 2.501945e+06 6.573359e+01 1.266710e+04

vif(ms)
         x1s          x2s          x3s        x1x2s        x1x3s        x2x3s         x12s 
2.856749e+06 1.095614e+04 2.017163e+06 9.802903e+03 1.428092e+06 2.403594e+02 2.501945e+06 
        x22s         x32s 
6.573359e+01 1.266710e+04 

---

Reference:

$\rm [I]$ Linear Regression Analysis, George A. F. Seber, Alan J. Lee, John Wiley & Sons, $2003,$ sec. $3.11.2,$ p. $71.$