R: Error in vif.default(m) : model contains fewer than 2 terms

Question

When I fit a polynomial regression model in R, if I use the function poly() and then try to get the variance inflation factors using vif() I get the following error:

y = c(0.22200,0.39500,0.42200,0.43700,0.42800,0.46700,0.44400,0.37800,0.49400,
      0.45600,0.45200,0.11200,0.43200,0.10100,0.23200,0.30600,0.09230,0.11600,
      0.07640,0.43900,0.09440,0.11700,0.07260,0.04120,0.25100,0.00002)
x1 = c(7.3,8.7,8.8,8.1,9.0,8.7,9.3,7.6,10.0,8.4,9.3,7.7,9.8,7.3,8.5,9.5,7.4,7.8,
       7.7,10.3,7.8,7.1,7.7,7.4,7.3,7.6)
x2 = c(0.0,0.0,0.7,4.0,0.5,1.5,2.1,5.1,0.0,3.7,3.6,2.8,4.2,2.5,2.0,2.5,2.8,2.8,
       3.0,1.7,3.3,3.9,4.3,6.0,2.0,7.8)
x3 = c(0.0,0.3,1.0,0.2,1.0,2.8,1.0,3.4,0.3,4.1,2.0,7.1,2.0,6.8,6.6,5.0,7.8,7.7,
       8.0,4.2,8.5,6.6,9.5,10.9,5.2,20.7)
m = lm(y~poly(x1, x2, x3, degree=2, raw=TRUE))
summary(m)

Now call vif()

> vif(m)
Error in vif.default(m) : model contains fewer than 2 terms

The model has 9 terms and an intercept.

> m$rank
[1] 10

It seems to me that the vif() function does not function with poly(). Is this correct? Is there a way around this or do I need to use first principles?

I can calculate the variance inflation factors like so:

X = poly(x1, x2, x3, degree=2, raw=TRUE)
C = solve(t(X)%*%X)
vifs = 1/diag(C)

Answer 1

I don't see much point in using poly with raw=TRUE, especially in this case where you want multiple terms and poly delivers a very opaque labeling of its results. It's also unclear which version of vif is being used. I chose to experiment with using I() to create the "pure" second degree terms ( ... is "homogenous" the correct term here?) and the R formula interface (using (...)^2) for the rest and had no trouble understanding the results in contrast to what I got with poly;

> m = lm(y~(x1+x2+x3)^2 + I(x1^2) + I(x2^2) + I(x3^2))
> summary(m)

Call:
lm(formula = y ~ (x1 + x2 + x3)^2 + I(x1^2) + I(x2^2) + I(x3^2))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.063213 -0.037282 -0.001113  0.016738  0.122539 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -1.769364   1.286976  -1.375   0.1881  
x1           0.420799   0.294173   1.430   0.1718  
x2           0.222453   0.130742   1.701   0.1082  
x3          -0.127995   0.070245  -1.822   0.0872 .
I(x1^2)     -0.019325   0.016797  -1.150   0.2668  
I(x2^2)     -0.007449   0.012048  -0.618   0.5451  
I(x3^2)      0.000824   0.001441   0.572   0.5754  
x1:x2       -0.019876   0.012037  -1.651   0.1182  
x1:x3        0.009151   0.007621   1.201   0.2473  
x2:x3        0.002576   0.007039   0.366   0.7192  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06092 on 16 degrees of freedom
Multiple R-squared:  0.9169,    Adjusted R-squared:  0.8702 
F-statistic: 19.63 on 9 and 16 DF,  p-value: 5.051e-07

> rms::vif(m)
       x1        x2        x3   I(x1^2)   I(x2^2)   I(x3^2)     x1:x2     x1:x3     x2:x3 
521.01297 401.58833 688.02220 501.50614 173.60055  99.67708 204.43081 456.00750 349.97018 

In an effort to figure out what version of vif was "not working", I tried both rms::vif and HH::vif and neither one threw the error you encountered so I'm not sure why you got the error:

> m = lm(y~poly(x1, x2, x3, degree=2, raw=TRUE))
> HH::vif(m)
poly(x1, x2, x3, degree = 2, raw = TRUE)1.0.0 poly(x1, x2, x3, degree = 2, raw = TRUE)2.0.0 
                                    521.01297                                     501.50614 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.1.0 poly(x1, x2, x3, degree = 2, raw = TRUE)1.1.0 
                                    401.58833                                     204.43081 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.2.0 poly(x1, x2, x3, degree = 2, raw = TRUE)0.0.1 
                                    173.60055                                     688.02220 
poly(x1, x2, x3, degree = 2, raw = TRUE)1.0.1 poly(x1, x2, x3, degree = 2, raw = TRUE)0.1.1 
                                    456.00750                                     349.97018 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.0.2 
                                     99.67708 
> rms::vif(m)
poly(x1, x2, x3, degree = 2, raw = TRUE)1.0.0 poly(x1, x2, x3, degree = 2, raw = TRUE)2.0.0 
                                    521.01297                                     501.50614 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.1.0 poly(x1, x2, x3, degree = 2, raw = TRUE)1.1.0 
                                    401.58833                                     204.43081 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.2.0 poly(x1, x2, x3, degree = 2, raw = TRUE)0.0.1 
                                    173.60055                                     688.02220 
poly(x1, x2, x3, degree = 2, raw = TRUE)1.0.1 poly(x1, x2, x3, degree = 2, raw = TRUE)0.1.1 
                                    456.00750                                     349.97018 
poly(x1, x2, x3, degree = 2, raw = TRUE)0.0.2 
                                     99.67708 

Maybe it was from car?

> car::vif(m)
Error in vif.default(m) : model contains fewer than 2 terms

If so the proper way to handle this would seem to be sending an issues submission to John Fox. Does he have a GitHub page for issues? You should be able to find the answer to that question by looking at the output packageDescription("car"). Nope. You will neeed to send him an email. Here's how to get that destination:

maintainer("car")

Addendum:

m = lm(y~(x1+x2+x3)^2 + I(x1^2) + I(x2^2) + I(x3^2))

> car::vif(m)
there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif
       x1        x2        x3   I(x1^2)   I(x2^2)   I(x3^2)     x1:x2     x1:x3     x2:x3 
521.01297 401.58833 688.02220 501.50614 173.60055  99.67708 204.43081 456.00750 349.97018 

So maybe there's some useful information in the help page? Might be worth seeing what that message is all about.