使用所有可能的var组合自动化lm测试，并获取R中的shapiro.test()、bptest()、vif()的值

Question

我花了几天的时间在R中寻找能够满足所有标准OLS假设(正态分布、同方差、无多重共线性)的最优模型，但在12个变量中，不可能找到最优的var组合。因此，我正在尝试创建一个脚本，它将自动执行此过程。

下面是用于计算的示例代码：

x1 <- runif(100, 0, 10)
x2 <- runif(100, 0, 10)
x3 <- runif(100, 0, 10)
x4 <- runif(100, 0, 10)
x5 <- runif(100, 0, 10)

df <- as.data.frame(cbind(x1,x2,x3,x4,x5))

library(lmtest)
library(car)

model <- lm(x1~x2+x3+x4+x5, data = df)

# check for normal distribution (Shapiro-Wilk-Test)
rs_sd <- rstandard(model)
shapiro.test(rs_sd)

# check for heteroskedasticity (Breusch-Pagan-Test)
bptest(model)

# check for multicollinearity
vif(model)

#-------------------------------------------------------------------------------
# models without outliers
# identify outliers (calculating the Cooks distance, if x > 4/(n-k-1) --> outlier
cooks <- round(cooks.distance(model), digits = 4)
df_no_out <- cbind(df, cooks)
df_no_out <- subset(df_no_out, cooks < 4/(100-4-1))

model_no_out <- lm(x1~x2+x3+x4+x5, data = df_no_out)

# check for normal distribution
rs_sd_no_out<- rstandard(model_no_out)
shapiro.test(rs_sd_no_out)

# check for heteroskedasticity
bptest(model_no_out)

# check for multicollinearity
vif(model_no_out)

我考虑的是遍历所有的var组合，并获得shapiro.test()和bptest()的P值或创建的所有模型的VIF值，以便我可以比较显着值或多重共线性。(在我的数据集中，多重共线性应该不是问题，因为为了检查多重共线性，VIF测试会产生更多的值(对于每个var 1xVIF因子)，这可能会对在代码中实现更具挑战性)，shapiro.test + bptest()的p值将足以满足…)。

我试着写了几个脚本来自动化这个过程，但是没有成功(不幸的是我不是程序员)。我知道已经有一些线程在处理这个问题。

How to run lm models using all possible combinations of several variables and a factor

Finding the best combination of variables for high R-squared values

但我还没有找到一个脚本，它也可以计算P值。

特别是，对于没有异常值的模型的测试是很重要的，因为在去除异常值之后，OLS假设在许多情况下都是完全满足的。

我真的非常感谢任何建议或帮助这一点。

Answer 1

你触及了现在所说的统计学习的皮毛。介绍文本是“统计学习在R中的应用”，研究生水平的文本是“统计学习的要素”。要执行所需的操作，请使用"leaps“包中的regsubsets()函数。然而，如果你至少阅读了简介书中的第6章，你就会发现交叉验证和自举是进行模型选择的现代方法。

Answer 2

下面将自动执行模型拟合和之后所做的测试。

有一个函数可以适用于所有可能的模型。然后，对*apply函数的一系列调用将获得您想要的值。

library(lmtest)
library(car)


fitAllModels <- function(data, resp, regr){
  f <- function(M){
    apply(M, 2, function(x){
      fmla <- paste(resp, paste(x, collapse = "+"), sep = "~")
      fmla <- as.formula(fmla)
      lm(fmla, data = data)
    })
  }
  regr <- names(data)[names(data) %in% regr]
  regr_list <- lapply(seq_along(regr), function(n) combn(regr, n))
  models_list <- lapply(regr_list, f)
  unlist(models_list, recursive = FALSE)
}

现在来看数据。

# Make up a data.frame to test the function above.
# Don't forget to set the RNG seed to make the
# results reproducible
set.seed(7646)
x1 <- runif(100, 0, 10)
x2 <- runif(100, 0, 10)
x3 <- runif(100, 0, 10)
x4 <- runif(100, 0, 10)
x5 <- runif(100, 0, 10)

df <- data.frame(x1, x2, x3, x4, x5)

首先，以"x1"作为响应，其他变量作为可能的回归变量，对所有模型进行拟合。可以使用一个响应和任意数量的可能回归函数来调用该函数。

fit_list <- fitAllModels(df, "x1", names(df)[-1])

现在是测试的顺序。

# Normality test, standardized residuals
rs_sd_list <- lapply(fit_list, rstandard)
sw_list <- lapply(rs_sd_list, shapiro.test)
sw_pvalues <- sapply(sw_list, '[[', 'p.value')

# check for heteroskedasticity (Breusch-Pagan-Test)
bp_list <- lapply(fit_list, bptest)
bp_pvalues <- sapply(bp_list, '[[', 'p.value')

# check for multicollinearity, 
# only models with 2 or more regressors
vif_values <- lapply(fit_list, function(fit){
  regr <- attr(terms(fit), "term.labels")
  if(length(regr) < 2) NA else vif(fit)
})

关于厨师距离的说明。在您的代码中，您正在设置原始data.frame的子集，生成一个没有异常值的新a。这将复制数据。我只选择了df行的索引列表。如果您更喜欢重复的data.frames，请取消注释下面匿名函数中的行，并注释掉最后一行。

# models without outliers
# identify outliers (calculating the 
# Cooks distance, if x > 4/(n - k - 1) --> outlier

df_no_out_list <- lapply(fit_list, function(fit){
  cooks <- cooks.distance(fit)
  regr <- attr(terms(fit), "term.labels")
  k <- length(regr)
  inx <- cooks < 4/(nrow(df) - k - 1)
  #df[inx, ]
  which(inx)
})

# This tells how many rows have the df's without outliers
sapply(df_no_out_list, NROW)

# A data.frame without outliers. This one is the one 
# for model number 8. 
# The two code lines could become a one-liner.
i <- df_no_out_list[[8]]
df[i, ]