Variance Inflation Factors (VIFs) on model vs covariates themselves

Question

I am confused on which type of "object" do the VIF functions operate.

Let me give two examples, which are confusing me. The VIFs from the car and AED libraries are purportedly doing a very similar thing. However:

1) The vif() command in the R package car calculates VIFs based on the model (for example, a linear model). I have no issue with interpreting the results.

But this is clearly different from the following:

2) Zuur et al. 2009 (Mixed effects models and extensions in ecology with R) have produced their corvif function within the AED package. There, the VIFs seem to be calculated based on the covariates themselves (i.e. before the model is even fitted).

Here is an example from their book (I am not including actual data here, but that's not the point anyway):

library(AED); data(Tbdeer)
Z <- cbind(Tbdeer$OpenLand, Tbdeer$ScrubLand,
Tbdeer$QuercusPlants, Tbdeer$QuercusTrees,
Tbdeer$ReedDeerIndex, Tbdeer$EstateSize,
Tbdeer$Fenced)
corvif(Z)

And that's what is confusing me. Also, the vif() command from the car package does not seem to work unless the object is a model.

(I realise this is somewhat related to R and coding, but I thought it has a more general statistical relevance, so I posted it here)

Any thoughts?

Answer 1

The Variance Inflation Factor (VIF) can be defined as

$$\frac{1}{1-R_i^2}$$ where $R_i^2$ is the R-squared value for the regression of the $i$-th regressor on the other regressors.

So all you need to know are the regressors and not the observed outcomes. The object on which the VIF functions operate is the regressor-matrix.

---

The R-code below shows four methods that all work the same and return a vector

15.373833 21.620241  9.832037  3.374620 15.164887 7.527958  4.965873  4.648487  5.357452  7.908747

The four methods are

1. The corvif function from the aed package (since this is discontinued a copy of the code is added below).

   This function uses just the regressor-matrix.

2. The vif function from the car package.

   This function also uses the regressor-matrix. But, it does this indirectly. You need to give an object to the function from which the regressor-matrix can be obtained using functions model.matrix, coef, or vcov.

3. The computation using $$\frac{1}{1-R_i^2}$$

4. An alternative computation with $$\text{Var}(X_j) \cdot [(X^tX)^{-1}]_{jj} \cdot (n-1)$$ where $n$ is the number of observations.

code:

### data and properties like number of parameters p and number of observations n
### 
data = datasets::mtcars
Y <- data$mpg
Z <- cbind(data$cyl, data$disp, data$hp, data$drat, data$wt, data$qsec, data$vs, data$am, data$gear, data$carb)
n = length(Y)
p = length(Z[1,])

### design matrix / regressor-matrix
Zp <- cbind(rep(1,n),Z)  
### linear model containing the design matrix
mod <- lm(Y ~ . , data = as.data.frame(Z))
          
### aed::corvif
corvif(Z)

### car::vif
car::vif(mod)

### computation with R-squared of regressing regressor vs other regressors
sapply(1:p, FUN = function(i) { 
  mod <- lm(Z[,i] ~ 1+Z[,-i])
  1/(1-summary(mod)$r.squared)
})

### computation with inverse of model matrix (X^tX)^-1
apply(Z,2,var) *diag(solve(t(Zp) %*% Zp))[-1] * (n-1)
      

---

Code for corvif

corvif <- function(dataz) {
  dataz <- as.data.frame(dataz)
  #correlation part
  #cat("Correlations of the variables\n\n")
  #tmp_cor <- cor(dataz,use="complete.obs")
  #print(tmp_cor)
  
  #vif part
  form    <- formula(paste("fooy ~ ",paste(strsplit(names(dataz)," "),collapse=" + ")))
  dataz   <- data.frame(fooy=1,dataz)
  lm_mod  <- lm(form,dataz)
  
  cat("\n\nVariance inflation factors\n\n")
  print(myvif(lm_mod))
}


#Support function for corvif. Will not be called by the user
myvif <- function(mod) {
  v <- vcov(mod)
  assign <- attributes(model.matrix(mod))$assign
  if (names(coefficients(mod)[1]) == "(Intercept)") {
    v <- v[-1, -1]
    assign <- assign[-1]
  } else warning("No intercept: vifs may not be sensible.")
  terms <- labels(terms(mod))
  n.terms <- length(terms)
  if (n.terms < 2) stop("The model contains fewer than 2 terms")
  if (length(assign) > dim(v)[1] ) {
    diag(tmp_cor)<-0
    if (any(tmp_cor==1.0)){
      return("Sample size is too small, 100% collinearity is present")
    } else {
      return("Sample size is too small")
    }
  }
  R <- cov2cor(v)
  detR <- det(R)
  result <- matrix(0, n.terms, 3)
  rownames(result) <- terms
  colnames(result) <- c("GVIF", "Df", "GVIF^(1/2Df)")
  for (term in 1:n.terms) {
    subs <- which(assign == term)
    result[term, 1] <- det(as.matrix(R[subs, subs])) * det(as.matrix(R[-subs, -subs])) / detR
    result[term, 2] <- length(subs)
  }
  if (all(result[, 2] == 1)) {
    result <- data.frame(GVIF=result[, 1])
  } else {
    result[, 3] <- result[, 1]^(1/(2 * result[, 2]))
  }
  invisible(result)
}
#END VIF FUNCTIONS