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,
x1 = c(1300,1300,1300,1300,1300,1300,1200,1200,1200,1200,1200,1200,
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,

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

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 

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.


2 Answers 2


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)


[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

         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 


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.