R:VIF 自定义函数

R: customised function for VIF

我正在尝试编写一个循环来计算方差 Inflation 因子。我知道有些功能和包可以为我做这件事,但我需要某种定制。

样本数据

  library(MASS)
  library(clusterGeneration)

  set.seed(2)
  num.vars <- 30
  num.obs<-200
  cov.mat<- genPositiveDefMat(num.vars,covMethod="unifcorrmat")$Sigma
  rand.vars<- mvrnorm(num.obs,rep(0,num.vars),Sigma=cov.mat)

  cov.mat <- as.data.frame(cov.mat)
  names(cov.mat) <- rep(paste0("X",1:30))

此数据框有 30 列(预测变量)。

这是我的循环逻辑:

1) 将每个预测变量与其他预测变量进行回归并计算 R2。使用 VIF = 1/1 - R2 将 R2 转换为 VIF。这将给我 30 个 VIF 值。

2) 对VIF值进行排序。如果顶部预测变量的 VIF > 10,则从 cov.mat 中删除预测变量。 cov.mat 现在将有 29 个预测变量。

3) 重复步骤 1,即将每个预测变量与其他预测变量进行回归并再次计算 VIF(这次是 29 个 VIF)。如果 max VIF > 10,删除具有最高 VIF 的变量并继续执行直到 max VIF <= 10.

但是,问题是我想保留 X4 、 X6 和 X10 ,即使它们在给定迭代中的 VIF > 10 也是如此。所以在上面的过程中,如果 X4 或 X6 或 X10 在一次迭代中得出具有最高 VIF(> 10),则删除具有第二高 VIF 的变量(仅当第二高 VIF 也 > 10 并且不是 X4 或X6 或 X10)。我希望这是清楚的

  mat <- matrix(, ncol = 2, nrow = nrow(cov.mat)) #  this will store the 30 VIFs

for(i in 1:ncol(cov.mat)){
      mdl <- lm(cov.mat[,i] ~ ., data = cov.mat) # this will regress each column against other columns but throws an error when i = 2
      r.squared <- unlist(summary(mdl)[8]) # this gives the r-squared of predictor i
      vif <- 1/(1- r.squared^2) # calcualtion of VIF for predictor i
      mat[i,2]  <- vif
      mat[i,1]  <- names(cov.mat[i])
  }

假设上面的循环工作正常,我有一个矩阵,第一列是变量名,第二列是 VIF 值。

     df <- data.frame(mat)
     names(df) <- c("variable", "vif")
     df <- df[sort(df$vif),]

     ifelse(df[1,2] <= 10, stop, ifelse(df[1,2] > 10 & names(df[1,1]) != "X4" | names(df[1,1]) != "X6" | names(df[1,1]) != "X10", ....

这就是我迷路的地方。

我首先需要检查具有最高 VIF 的变量是否 > 10 并且不在 X4 或 x6 和 X10 之间,然后从数据帧 cov.mat 中删除变量。 如果具有最高 VIF 的变量(给定 VIF > 10)是 X4 或 X6 或 X10,则转到 df 的第二行并评估其 VIF 是否 > 10 和 是否不在X4、X6、X10中,如果满足则从cov.mat中移除,重新开始迭代。

编辑

我的原始数据框有 51 列和 1458 行。当我运行上面的函数时,它给我一个错误there are aliased coefficients in the model。 为什么会这样?

在您的示例数据中,无法计算整个数据集的 或 VIF 分数,很可能是因为完全共线性。但是,此处的函数应该适用于不是这种情况的数据(例如,数据集的列 1:15)。您可以 ignore/remove 所有 cat 代码。那只是为了说明正在发生的事情

此外,我使用包 car 作为函数 vif

library(vif)

vif_fun <- function(df, keep_in) {
             # df: the dataset of interest
             # keep_in: the variables that should be kept in  
             highest <- c()
             while(TRUE) {
                # the rnorm() below is arbitrary as the VIF should not 
                # depend on it
                vifs <- vif(lm(rnorm(nrow(df)) ~. , data = df))
                adj_vifs <- vifs[-which(names(vifs) %in% keep_in)]
                if (max(adj_vifs) < 10) {
                     break
                }
               cat("\n")
               print(vifs)
               highest <- c(highest,names((which(adj_vifs == max(adj_vifs)))))
               cat("\n")
               cat("removed:", highest)
               cat("\n")
               df <- df[,-which(names(df) %in% highest)]

              }
            cat("\n")
            cat("final variables: \n")
            return(names(vifs))
              }

# example with mtcars dataset
vif_fun(mtcars,keep_in = c("cyl"))


# example using part of your data
vif_fun(cov.mat[,1:15], keep_in = c("X15", "X12"))