代码之家 › 专栏 › 技术社区 › 89_Simple

R: VIF定制功能

r-car correlation regression r

89_Simple · 技术社区 · 6 年前

我试图编写一个循环来计算方差膨胀系数。我知道有一些功能和软件包可以为我做到这一点,但我需要一些定制。

A样本数据

  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。小地毯 现在将有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之间,并从数据帧中删除该变量 cov。小地毯 . 如果VIF最高的变量(给定VIF>10)为X4、X6或X10,则转到 df 并评估其VIF是否;10或10,以及是否不在X4、X6或X10之间,如果满足条件,则将其从 cov。小地毯 然后再次开始迭代。

编辑

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

1 回复 | 直到 6 年前

Daniel 6 年前

在您的示例数据中,无法计算整个数据集的or 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"))