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

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.

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

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