0

我花了几天时间寻找能够满足 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-VALUES 或创建的所有模型的 VIF 值，以便我可以比较显着性值或多重共线性 （在我的数据集中，多重共线性应该不是问题，因为要检查多重共线性，VIF 测试会产生更多值（对于每个 var 1xVIF 因子），这对于在代码中实现可能更具挑战性），p 值shapiro.test + bptest() 就足够了……）。

我尝试编写几个脚本来自动化该过程但没有成功（不幸的是我不是程序员）。我知道已经有一些线程在处理这个问题

如何使用多个变量和一个因子的所有可能组合运行 lm 模型

为高 R 平方值寻找最佳变量组合

但我还没有找到一个也可以计算 P-VALUES 的脚本。

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

我非常感谢任何建议或帮助。

4

2 回答 2

2

您正在触及现在称为统计学习的表面。介绍文本是“在 R 中应用的统计学习”，研究生级别的文本是“统计学习的要素”。要执行您需要的操作，请使用“跳跃”包中的 regsubsets() 函数。但是，如果您至少阅读了介绍书中的第 6 章，您会发现交叉验证和引导程序是进行模型选择的现代方式。

于 2018-11-16T15:06:25.027 回答
1

以下自动化模型拟合和您之后进行的测试。

有一个函数适合所有可能的模型。然后对函数的一系列调用*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 进行子集化，生成一个没有异常值的新数据。这将重复数据。我只选择了 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, ]
于 2018-11-16T16:26:49.353 回答