3
$\begingroup$

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?

$\endgroup$
9
  • $\begingroup$ This the explanation in the website: "In the book we use the package AED to load data. However, we haven given up compiling a new version of the AED package each time a new R version comes out. Therefore we no longer provide AED." Is there any other way to get the package? $\endgroup$
    – T.E.G.
    Feb 15, 2017 at 22:23
  • $\begingroup$ Yes, on the same page (highstat.com/book2.htm), they provide the code, saying: "To run the corvif function and the pairplot with the Pearson correlation coefficients, download the file HighstatLibV6.R (use right-mouse click and Save As), save it to your computer and in R type: > source("C:/YourDirectory/HighstatLibV6.R")." highstat.com/Book2/HighstatLibV6.R $\endgroup$
    – Tilen
    Feb 15, 2017 at 22:28
  • 2
    $\begingroup$ The functions return same results. You don't need to fit the model to calculate vif, if that was your question. You only need the list of independent variables. $\endgroup$
    – T.E.G.
    Feb 15, 2017 at 22:42
  • $\begingroup$ Yes, that was my question. But when I tried that via "vif" function in "car" library, I got an error. It says "Error: $ operator is invalid for atomic vectors". $\endgroup$
    – Tilen
    Feb 15, 2017 at 22:59
  • $\begingroup$ I think that is the coding part of the question. The documentation of the vif in car shows, "an object that responds to coef, vcov, and model.matrix, such as an lm or glm object" should be used. That might not be the case for other similar/equivalent functions. It is the authors' decision about how to write the functions. $\endgroup$
    – T.E.G.
    Feb 15, 2017 at 23:03

1 Answer 1

2
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.