Quadratic term and variance inflation factor in OLS estimation

Question

One question I have carried around with me for a while is related to including quadratic terms for model specifications. I wonder why it is considered ok to include linear and quadratic terms into OLS analysis, despite the fact that the Variance Inflation Factor (VIF) gets high. A model specification that I was taught to be fine to use is the following:

$y = c + x + x^2 + \epsilon$

But this model has quite high numerical linear correlation between $x$ and $x^2$ and I do not understand why in this case the VIF does not matter.

An example for what I mean using R is the following:

# Load packages
require(car)
requre(lmtest)

# Create artificial data
dat <- data.frame(x = 1:1000, x2 = (1:1000)^2)
dat$y <- dat$x2 + rnorm(1000, 0, 10)

# Run model with linear term
dat.lm1 <- lm(y ~ x, dat)

# Run model with linear and quadratic term
dat.lm2 <- lm(y ~ x + x2, dat)

In this setup, the RESET test suggests using quadratic terms but the VIF is high.

> reset(dat.lm1)
RESET test
data:  dat.lm1 
RESET = 27344126617, df1 = 2, df2 = 996, p-value < 2.2e-16

> reset(dat.lm2)
RESET test
data:  dat.lm2 
RESET = 0.4404, df1 = 2, df2 = 995, p-value = 0.6439

> vif(dat.lm2)
x       x2 
16.03008 16.03008 

So the question is: Is this really a model specification that is valid and if yes, why does the VIF not matter in that case?

Answer 1

Generally collinearity only affects variance estimates of your parameters. You can avoid collinearity by data transformation. For your example:

# norm x to mean zero 
dat$x_n <- dat$x - mean(dat$x)
dat$x2_n <- (dat$x_n)^2
dat.lm2_n <- lm(y ~ x_n + x2_n, dat)

This transformation doesn't affect your residuals but you get smaller confidence intervals.

> summary(dat.lm2)

Call:
lm(formula = y ~ x + x2, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-37.505  -7.182  -0.146   7.016  36.599 

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  5.967e-01  9.626e-01      0.620    0.535    
x           -1.185e-03  4.441e-03     -0.267    0.790    
x2           1.000e+00  4.296e-06 232752.480   <2e-16 ***
---
Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1 

Residual standard error: 10.13 on 997 degrees of freedom
Multiple R-squared:     1,  Adjusted R-squared:     1 
F-statistic: 4.342e+11 on 2 and 997 DF,  p-value: < 2.2e-16 

> summary(dat.lm2_n)

Call:
lm(formula = y ~ x_n + x2_n, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-37.505  -7.182  -0.146   7.016  36.599 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.505e+05  4.804e-01  521493   <2e-16 ***
x_n         1.001e+03  1.109e-03  902349   <2e-16 ***
x2_n        1.000e+00  4.296e-06  232752   <2e-16 ***
---
Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1 

Residual standard error: 10.13 on 997 degrees of freedom
Multiple R-squared:     1,  Adjusted R-squared:     1 
F-statistic: 4.342e+11 on 2 and 997 DF,  p-value: < 2.2e-16 

And the VIF shrinks of course:

> vif(dat.lm2)
       x       x2 
16.03008 16.03008 

> vif(dat.lm2_n)
 x_n x2_n 
   1    1