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"))
我正在尝试编写一个循环来计算方差 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"))