R: индивидуальная функция для VIF

Question

Я пытаюсь написать цикл для расчета фактора инфляции дисперсии. Я понимаю, что есть функции и пакеты, которые могут сделать это для меня, но мне нужна какая-то настройка.

Образец данных

  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. Конвертировать R2 в VIF, используя VIF = 1/1 - R2. Это даст мне 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 или Х6 или Х10). Надеюсь это понятно

  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, Почему это происходит?

Answer 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"))