Multicollinearity and VIF: a worked example in R

Multicollinearity is overlap between the predictors in a regression: two or more of them carry much of the same information. It does not bias the coefficients, it makes them imprecise, and the variance inflation factor (VIF) says by how much. Whether that matters depends on what the model is for.

The post on linear regression assumptions covers collinearity in a few paragraphs among the other checks. This one is the dedicated worked example: a real regression on R’s swiss data, where the VIF is computed from first principles and the overlap is moderate, then a simulated pair of predictors where it is extreme, so you can see exactly what a high VIF does to a confidence interval. Every example number is a deterministic engine output or arithmetic on one.

What multicollinearity is, and what it is not

In a multiple regression each coefficient is the change in the mean outcome for a one-unit step in that predictor, holding the other predictors fixed. That “holding fixed” is the trouble. If two predictors rise and fall together, the data contain very few cases where one differs and the other does not, so the model has little information to separate their effects.

The consequences are narrow and specific:

• The estimates stay unbiased. In a correctly specified model, overlap does not push the coefficients in any direction on average.
• The standard errors grow. The intervals for the overlapping coefficients widen, and a coefficient can look unimportant even when the predictors jointly matter a lot.
• The estimates become unstable. Across samples, or when a correlated predictor is added or removed, an individual coefficient can move a long way, and can change sign.
• The model as a whole is not harmed. The overall fit, and the predictions for cases that follow the same pattern of overlap, are not degraded.

Two terms are worth separating. Collinearity is overlap between two predictors; multicollinearity is the same idea for several, where one predictor is nearly a combination of a few others even if no single pair is strongly correlated. A correlation matrix finds the first and can miss the second, which is why a per-predictor measure is used.

At the extreme, perfect collinearity, one predictor is an exact combination of others (percentages that sum to 100, or a variable entered twice in different units). Then the coefficients cannot be separated at all, and R’s lm() returns NA for the redundant coefficient (see the lm documentation). That is a data problem to fix, not a statistical one to manage.

How the VIF is computed

For each predictor, regress it on all the other predictors, ignoring the outcome. The R-squared of that auxiliary regression says how much of the predictor the others already explain. The VIF is

VIF = 1 / (1 - R-squared)

so a predictor that the others do not explain at all has a VIF of 1, and the VIF climbs without limit as the R-squared approaches 1. It is the factor by which the variance of that coefficient is larger than it would be if the predictor were uncorrelated with the rest, and its square root is the same factor for the standard error and the width of the confidence interval (O’Brien 2007).

Three consequences follow from the definition:

1. The VIF is a property of the predictors alone. The outcome does not enter, so it does not change if you swap the outcome.
2. It is one number per coefficient, not per model. Some coefficients can be affected while others are untouched.
3. It measures inflation relative to a no-overlap baseline. It does not say whether the inflated interval is too wide for your purpose.

For a factor predictor with several levels, one coefficient per level makes a single VIF meaningless. The generalized VIF handles the whole set of terms together (Fox and Monette 1992), and GVIF^(1/(2*Df)) puts it on the same scale as the square root of an ordinary VIF, so the same comparisons apply.

What counts as a high VIF

You will see 5 and 10 quoted as limits, and sometimes 4. They are conventions. O’Brien (2007) examined them and concluded that a VIF above such a threshold is not by itself a reason to question a result or to drop a predictor, and that the usual remedies can do more harm than the condition they treat. The practical question is whether the confidence interval for the coefficient you care about is narrow enough for the decision in front of you.

A VIF of 10 means the standard error is about 3.2 times larger than with no overlap; a VIF of 4 means 2 times. Read it as a multiplier on interval width, and judge the width. Belsley, Kuh and Welsch (1980) is the standard book-length treatment of collinearity diagnostics beyond the VIF.

Worked example: fertility in 47 Swiss provinces

R’s swiss data hold one row for each of 47 French-speaking Swiss provinces around 1888, with a standardised fertility measure and five socio-economic variables: the share of males working in agriculture (Agriculture), the share of army draftees with the top mark on an exam (Examination), the share with education beyond primary school (Education), the share Catholic (Catholic) and infant mortality (Infant.Mortality). The research question for the run: how is fertility related to these variables across provinces? The rows are treated as independent provinces and the outcome is continuous, so a multiple linear regression fits. These are province-level associations in observational data, not statements about individuals or about what would happen if a province changed.

The engine fitted the model on all 47 rows (41 residual degrees of freedom). The predictors together give R-squared = 0.7067.

Predictor	Estimate	95% CI	p
Agriculture	-0.1721	-0.3141 to -0.03013	0.0187
Examination	-0.2580	-0.7707 to 0.2547	0.315
Education	-0.8709	-1.2406 to -0.5013	< 0.001
Catholic	0.1041	0.03291 to 0.1753	0.0052
Infant.Mortality	1.077	0.3061 to 1.848	0.0073

The table holds five coefficient tests with no multiplicity correction, so each p-value describes its coefficient on its own and the borderline ones should not be over-read.

The Examination row is the one to look at. Its interval runs from -0.7707 to 0.2547, wide enough to be compatible with a sizeable negative association and with a small positive one. The next step is to find out how much of that width is overlap.

Computing the VIF for Examination

A second engine run regressed Examination on the other four predictors. Its R-squared was 0.7279: nearly three quarters of the variation in Examination across provinces is explained by Agriculture, Education, Catholic and Infant.Mortality, with Education and Catholic the two clearly associated terms in that auxiliary model (both p < 0.001). So

VIF = 1 / (1 - 0.7279) = 3.68

and the standard error of the Examination coefficient is about 1.9 times what it would be with no overlap (the square root of 3.68). The engine’s own screen, which flags any predictor above a VIF of 5, raised no collinearity note for this model, so none of the five VIFs exceeded 5.

That is a typical outcome for observational data: a VIF under 5 and still an interval nearly twice as wide as it could have been. In R the same quantities come from one call:

fit <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality,
          data = swiss)
car::vif(fit)

# the same number for one predictor, by hand
aux <- lm(Examination ~ Agriculture + Education + Catholic + Infant.Mortality, data = swiss)
1 / (1 - summary(aux)$r.squared)

What happens if you drop the overlapping predictor

The common advice is to remove one of the overlapping variables. Here is what that does. A third run refitted the model without Education:

Predictor	Estimate	95% CI	p
Agriculture	-0.04759	-0.2097 to 0.1145	0.557
Examination	-0.9680	-1.4783 to -0.4578	< 0.001
Catholic	0.02611	-0.05144 to 0.1037	0.501
Infant.Mortality	1.396	0.4624 to 2.330	0.0043

Three things changed, and none of them is “the problem went away”.

• The Examination estimate moved from -0.2580 to -0.9680. The coefficient now carries some of what Education carried. The model answers a different question: the association of examination results with fertility, not holding education fixed.
• The interval did not get narrower. It is about as wide as before (1.02 against 1.03 units), because removing a predictor also raises the residual variation, which offsets the smaller inflation. What changed is where it sits: it now excludes zero.
• The model explains less. R-squared fell from 0.7067 to 0.5448, and the other coefficients moved too, Agriculture and Catholic becoming indistinguishable from zero.

One caution on this refit: its residuals show evidence of non-normality (Shapiro-Wilk p = 0.022, where the full model’s residuals showed none, p = 0.932), so its p-values and intervals deserve less trust than the full model’s. Neither version is “the right one” by this comparison alone. If the question was about examination results while taking education into account, the full model is the one that asks it, with an honest, wide interval. Dropping Education to make Examination look significant would be choosing the answer by choosing the model.

A high VIF: two predictors that correlate at 0.96

The swiss overlap is moderate. To see a high VIF, the engine was run on simulated data (100 independent rows, seeded, nothing real). The outcome was built as y = 1 + 0.5 x1 + 0.5 x2 + error with unit error spread, so both predictors truly matter and the true coefficients are known. In one dataset x2 was built to correlate with x1 at about 0.95 (realised: 0.96); in the other x2 was generated independently. Everything else, including the seed and the error draws, was the same.

The correlated fit:

Predictor	Estimate	95% CI	p
x1	1.013	0.3927 to 1.633	0.0016
x2	-0.01934	-0.6355 to 0.5969	0.950

The engine flagged both predictors for a VIF above 5. The auxiliary regression of x1 on x2 gave R-squared = 0.9250, so

VIF = 1 / (1 - 0.9250) = 13.3

and the standard errors should be about 3.65 times those of a no-overlap design. The independent version is the comparison:

Predictor	Estimate	95% CI	p
x1	0.5194	0.3493 to 0.6896	< 0.001
x2	0.3378	0.1454 to 0.5302	0.0007

The standard error of x1 was 0.3124 in the correlated data and 0.08572 in the independent data, a ratio of 3.6, as the VIF predicted. The x1 interval is 1.24 units wide against 0.34. Both intervals contain the true value of 0.5, so nothing is biased; the correlated data just cannot pin the value down. The estimate for x2 is about zero with an interval from -0.6355 to 0.5969, even though its true effect equals that of x1. Read alone, that row would say “no evidence x2 matters”. The correct reading is that the data cannot separate x1 from x2.

This is the situation in which a high VIF is a real problem for inference about individual coefficients. The interval widths, and the unstable and even misleading point estimates, are the damage the number describes.

When multicollinearity is not a problem

A high VIF only hurts the coefficients it is attached to, and only for the questions that need them.

• The collinear variables are controls, not the question. The VIF of your main predictor depends on its own overlap with the others. If two adjustment variables are tightly correlated with each other but not with the predictor you care about, the interval for that predictor is unaffected by their overlap.
• The question is about the group, not the individuals. If the overlapping predictors matter jointly, test them together: a joint F test comparing the model with and without the whole set does not suffer from their overlap, and in R it is anova(reduced_fit, full_fit).
• You want predictions. Predictions for new cases with the same pattern of overlap are not degraded by it. Predictions for cases that break the pattern extrapolate and are unreliable.
• The overlap is built in. A squared term and its linear term, or an interaction and its components, are correlated by construction. Centering the predictors changes what the lower-order coefficients mean, but it does not change the model’s fit or the precision of the product term’s estimate (Echambadi and Hess 2007), so a high VIF there is not evidence that anything is wrong.

What to do about a high VIF that matters

When inference about individual overlapping coefficients is the goal, here are the options, in order of how little they change the question:

1. Report the interval as it is. A wide interval is an honest statement of what the data can say. Reporting the VIF next to it tells the reader why.
2. Test the overlapping predictors together when the question allows it, as above.
3. Combine predictors that measure one thing into a single index or scale, decided on subject-matter grounds before you look at the outcome. Two exam-related variables, or two measures of the same exposure, are candidates; a combination chosen to make the result significant is not.
4. Get more data, or data in which the predictors vary independently. The variance of a coefficient depends on the sample size as well as the VIF, so a larger sample offsets the inflation (O’Brien 2007).
5. Drop a predictor only when it is redundant for the question, and say so, as the swiss comparison showed: the model you fit afterwards is a different model.

Shrinkage methods such as ridge regression, and principal components, are other routes. They trade interpretability or some bias for stability, so they answer a modified question and are best chosen deliberately, not as a way to clear a diagnostic (O’Brien 2007).

How to report it

State the model, the VIFs or the largest one, and the consequence for the coefficients that matter:

Variance inflation factors were examined for each predictor. In the swiss regression the Examination coefficient had a VIF of 3.68 (auxiliary R-squared = 0.7279), so its standard error was about 1.9 times what it would be with uncorrelated predictors; its interval, -0.7707 to 0.2547, is correspondingly wide.

Do not write that multicollinearity “was not a problem” because every VIF was under 10, or “was a problem” because one was over 5. Quote the number and the width of the interval, and let the reader judge.

Choosing and running the example

This is the multiple linear regression covered in the post on multiple linear regression and when to switch, and how to choose a statistical test starts from the outcome type if you are unsure that a linear model fits. The linear regression method reference lists what the engine computes for it. Confidence intervals are read as in the overview of p-values and confidence intervals: the range of values the data are compatible with, not a probability about the coefficient.

The example’s reproduce.R reads the bundled swiss data and refits the five-predictor model with lm(), so it recomputes the estimates, intervals, p-values and R-squared in the first table, but not the VIF, the other runs, or the engine’s checks. The car::vif() call and the auxiliary regression above compute the VIF from the same data.

# Results

- Test: Linear regression
- n: 47
- statistic: -2.448
- p-value: 0.0187
- effect size: R-squared = 0.7067
- 95% CI: [-0.3141, -0.03013]
- seed: 20261023

## Assumption checks
- normality (Shapiro-Wilk): passed=TRUE (no evidence against normality)

## Power (sensitivity, not observed power)
- minimum detectable effect at n=47: Cohen's f2 = 0.3118
- to detect a conventional medium effect: 92 total to detect f2=0.15
- alpha = 0.05, target power = 0.8

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Linear regression  |  engine 0.31.1  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261023)
.data_path <- "data/multicollinearity-vif.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("Fertility", "Agriculture", "Examination", "Education", "Catholic", "Infant.Mortality")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- summary(lm(Fertility ~ Agriculture + Examination + Education + Catholic +      Infant.Mortality, data = data))
print(result)

# Sensitivity: minimum detectable effect at n=47 (power=0.8, alpha=0.05). NOT observed power.
power <- pwr::pwr.f2.test(u = 5, v = 41, sig.level = 0.05, power = 0.8)
print(power)