1

我正在尝试编写一个循环来计算方差通货膨胀因子。我知道有一些功能和包可以为我做到这一点，但我需要某种定制。

样本数据

  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）。如果最大 VIF > 10，则删除具有最高 VIF 的变量并继续执行直到最大 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 之间，然后从 dataframe 中删除该变量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。为什么会这样？

4

1 回答 1

1

在您的示例数据中，无法为整个数据集计算 或 VIF 分数，这很可能是因为完美的共线性。然而，这里的函数应该适用于不是这种情况的数据（例如，数据集的第 1:15 列）。您可以忽略/删除所有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"))
于 2018-01-23T21:52:59.073 回答