I kept running into the variance inflation factor used as a pass/fail gate: fit the model, call vif(), drop anything above 10. I wanted to know what that number actually buys you, so I put it on real data and measured. The short answer is that VIF tells you exactly one thing, the factor by which collinearity widens a standard error, and it cannot tell you on its own whether that widening matters. Two of the largest VIFs in this post turn out not to be collinearity at all.
The data is the EPA’s fuel economy database, a keyless zipped CSV covering every vehicle configuration the agency has rated. I keep gasoline vehicles from model year 2015 onward, one row per distinct configuration.
library(tidyverse)
library(car)
zip <- tempfile(fileext = ".zip")
download.file("https://www.fueleconomy.gov/feg/epadata/vehicles.csv.zip", zip, quiet = TRUE)
vehicles <- read_csv(unzip(zip, "vehicles.csv", exdir = tempdir()), show_col_types = FALSE)
cars <- vehicles |>
filter(year >= 2015, fuelType1 %in% c("Regular Gasoline", "Premium Gasoline"),
displ > 0, cylinders > 0, comb08 > 0, !is.na(drive)) |>
distinct(make, model, year, displ, cylinders, drive, trany, .keep_all = TRUE) |>
mutate(drive = fct_lump_min(factor(drive), 200))
nrow(cars)
## [1] 13115
That leaves 13,115 cars spanning model years 2015 to 2027. Engine displacement and cylinder count are the textbook collinear pair, and they do not disappoint: their correlation here is 0.92.
What does a high VIF actually cost?
A VIF of $v$ multiplies that coefficient’s standard error by $\sqrt{v}$, and does nothing else. It does not bias the estimate, shrink $R^2$, or damage predictions. The full standard error is
$$\text{SE}(\hat\beta_j) = \frac{\sigma}{s_j\sqrt{n-1}} \times \sqrt{\text{VIF}_j}$$
so VIF is one of four ingredients, alongside the residual spread $\sigma$, the predictor’s own spread $s_j$, and the sample size $n$. VIF knows nothing about the other three, which is exactly why it cannot tell you whether you have a problem.
m_both <- lm(comb08 ~ displ + cylinders, data = cars)
vif(m_both)
## displ cylinders
## 6.702701 6.702701
# the identity, checked against lm's own standard error
sigma(m_both) / (sd(cars$displ) * sqrt(nrow(cars) - 1)) * sqrt(vif(m_both)[["displ"]])
## [1] 0.08105514
coef(summary(m_both))["displ", "Std. Error"]
## [1] 0.08105514
Both terms come back at VIF 6.7, so their standard errors are 2.59 times wider than they would be if displacement and cylinder count were uncorrelated. That is the entire cost. If you want one block that runs on its own, the same calculation on a built-in dataset:
library(car)
vif(lm(mpg ~ disp + cyl + hp, data = mtcars))
## disp cyl hp
## 5.521460 6.732984 3.350964
Is a VIF of 6.7 a problem?
You cannot tell from the VIF. To show why, I refit the same model on 2,000 random subsamples at each of five sample sizes. Collinearity is a property of the predictors, so the VIF barely moves between them; only $n$ changes.
set.seed(42)
sizes <- c(60, 120, 250, 600, 2000)
sim <- map_dfr(sizes, \(n) map_dfr(1:2000, \(i) {
fit <- lm(comb08 ~ displ + cylinders, data = slice_sample(cars, n = n))
co <- coef(summary(fit))
tibble(n = n, vif = vif(fit)[["cylinders"]],
est = co["cylinders", "Estimate"], p = co["cylinders", "Pr(>|t|)"])
}))
summary_tbl <- sim |>
group_by(n) |>
summarise(median_vif = median(vif), sd_estimate = sd(est),
pct_wrong_sign = 100 * mean(est > 0), pct_significant = 100 * mean(p < 0.05))
summary_tbl
## # A tibble: 5 × 5
## n median_vif sd_estimate pct_wrong_sign pct_significant
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 60 6.99 0.821 12.2 15.7
## 2 120 6.92 0.566 5.3 28.5
## 3 250 6.78 0.366 1 52.4
## 4 600 6.74 0.237 0.05 90.1
## 5 2000 6.70 0.117 0 100
The median VIF sits between 6.7 and 7 at every sample size, and the consequences are nothing alike. At 60 cars the cylinder coefficient comes out with the wrong sign in 12.2% of samples and reaches significance in only 15.7%; at 2,000 cars it is never wrong and always significant. An identical VIF is harmless at one sample size and ruinous at another, which is the case against reading it as a threshold.
dsp_colors <- c("#0066CC", "#E8862D", "#159A6C", "#7D5BD6",
"#D64580", "#2AA9B8", "#C9A227")
dsp_theme <- theme_minimal(base_size = 13) +
theme(plot.background = element_rect(fill = "#ECECEF", color = NA),
panel.background = element_rect(fill = "#ECECEF", color = NA),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_line(color = "grey78"),
axis.ticks = element_blank(),
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold"))
band <- sim |>
group_by(n) |>
summarise(med = median(est), lo = quantile(est, .025), hi = quantile(est, .975),
wrong = 100 * mean(est > 0))
ggplot(band, aes(n, med)) +
geom_hline(yintercept = 0, linewidth = .6, color = "grey40") +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = dsp_colors[1], alpha = .22) +
geom_line(color = dsp_colors[1], linewidth = 1) +
geom_point(color = dsp_colors[1], size = 2.4) +
geom_text(aes(y = hi, label = paste0(sprintf("%.1f", wrong), "% wrong sign")),
vjust = -0.9, size = 3.4, color = "grey30",
data = filter(band, n <= 120)) +
scale_x_log10(breaks = sizes, labels = scales::comma,
expand = expansion(mult = .12)) +
scale_y_continuous(expand = expansion(mult = c(.05, .18))) +
labs(title = "Same VIF, very different cost",
subtitle = "Cylinder coefficient over 2,000 subsamples per size; VIF is ~6.7 throughout",
x = "Cars in the sample (log scale)", y = "mpg per extra cylinder") +
dsp_theme
The estimate is centered in the right place at every size. That is the substantive result worth keeping: holding engine displacement fixed, each additional cylinder costs about 0.81 mpg, and a study of 60 cars would report that effect backwards roughly one time in 8.
Should I drop one of the collinear predictors?
Dropping one buys almost nothing in prediction, and it quietly changes the question the model answers.
m_drop <- lm(comb08 ~ displ, data = cars)
c(both = coef(m_both)[["displ"]], dropped = coef(m_drop)[["displ"]])
## both dropped
## -2.390812 -3.477673
set.seed(7)
rmse <- map_dfr(1:200, \(i) {
idx <- sample(nrow(cars), floor(.7 * nrow(cars)))
train <- cars[idx, ]; test <- cars[-idx, ]
tibble(both = sqrt(mean((test$comb08 - predict(lm(comb08 ~ displ + cylinders, train), test))^2)),
dropped = sqrt(mean((test$comb08 - predict(lm(comb08 ~ displ, train), test))^2)))
})
colMeans(rmse)
## both dropped
## 4.586671 4.623335
Out-of-sample RMSE over 200 train/test splits is 4.59 mpg with both terms and 4.62 without, a difference of 0.04 mpg. But the displacement coefficient moves from -2.39 to -3.48, a 45% shift, because it stops being the effect of a liter at fixed cylinder count and becomes the effect of a liter with cylinder count free to rise alongside it. Both are correct answers to different questions, and dropping a predictor to satisfy a VIF rule silently swaps one for the other.
Why did my VIF explode after I added an interaction?
Because a raw interaction or squared term is a scaling artifact, not collinearity, and centering removes it without changing the model at all.
m_raw <- lm(comb08 ~ displ * drive, data = cars)
vif(m_raw)
## GVIF Df GVIF^(1/(2*Df))
## displ 7.972565 1 2.823573
## drive 8299.778987 4 3.089465
## displ:drive 10989.195481 4 3.199785
cars_c <- mutate(cars, displ_c = displ - mean(displ))
m_cen <- lm(comb08 ~ displ_c * drive, data = cars_c)
vif(m_cen)
## GVIF Df GVIF^(1/(2*Df))
## displ_c 7.972565 1 2.823573
## drive 3.743619 4 1.179401
## displ_c:drive 18.686738 4 1.441922
max(abs(fitted(m_raw) - fitted(m_cen)))
## [1] 4.650857e-11
The raw model reports a GVIF of 8,300 for drive and 11,000 for the interaction. Centering displacement drops those to 3.7 and 18.7. The two models produce fitted values that agree to 5e-11, the same $R^2$ and the same predictions: nothing about the fit changed, only the origin of the displacement axis. A squared term behaves the same way, with VIF falling from 23.1 to 1.7 on centering, and the reason is visible if you just plot the two terms against each other.
terms <- bind_rows(
tibble(x = cars$displ, y = cars$displ^2, panel = sprintf("Raw: r = %.3f", cor(cars$displ, cars$displ^2))),
tibble(x = cars_c$displ_c, y = cars_c$displ_c^2, panel = sprintf("Centered: r = %.3f", cor(cars_c$displ_c, cars_c$displ_c^2)))
) |>
mutate(panel = fct_inorder(panel))
ggplot(terms, aes(x, y)) +
geom_point(alpha = .10, color = dsp_colors[1], size = .8) +
facet_wrap(~panel, scales = "free") +
labs(title = "Centering is what kills the squared-term VIF",
subtitle = "Same model, same predictions: only the correlation between the two terms changes",
x = "Linear term", y = "Squared term") +
dsp_theme
Over the raw positive range the squared term is nearly a straight function of the linear one, so they correlate at 0.978. Centering makes it a U and the correlation falls to 0.648. As of car 3.1.5, vif() also prints a message on models like this suggesting type = "predictor", which aggregates each predictor with the terms it interacts with and returns a VIF of exactly 1 for both displ and drive here. That is the cleaner answer, and it says the same thing: the interaction never was collinear.
How do I read vif() when the model has a factor?
For any term with more than one degree of freedom, vif() returns a generalized VIF, and the column you compare is the third one, GVIF^(1/(2*Df)).
vif(lm(comb08 ~ displ + cylinders + drive, data = cars))
## GVIF Df GVIF^(1/(2*Df))
## displ 7.202009 1 2.683656
## cylinders 6.893096 1 2.625471
## drive 1.334533 4 1.036731
That last column is on the scale of a standard error multiplier rather than a variance multiplier, so the threshold matching the familiar VIF of 10 is $\sqrt{10} = 3.16$, not 10. The gap matters: in the interaction model above, drive shows a GVIF of 8,300 across 4 degrees of freedom, which reduces to 3.09, just under the conventional line. Reading the first column against 10 condemns a variable that the correct column says is fine.
None of this makes collinearity harmless. It makes VIF a narrower instrument than its reputation suggests: a standard error multiplier, meaningful only next to your sample size, and misleading the moment your model contains a squared term, an interaction, or a factor.
