在线性回归中使用数据框的列名作为预测变量

Question

我正在尝试遍历我的 data.frame 的所有列名并使用它们 作为线性回归中的预测变量。

我目前拥有的是：

for (i in 1:11){
for (j in 1:11){
if (i != j ){
  var1 = names(newData)[i]
  var2 = names(newData)[j]
  glm.fit = glm(re78 ~  as.name(var1):as.name(var2), data=newData)
  summary(glm.fit)
  cv.glm(newData, glm.fit, K = 10)$delta[1]
  }
 }
}

newData 是我的 data.frame，总共有 11 列。这段代码给了我以下错误：

model.frame.default(formula = re78 ~ as.name(var1), data = newData, 中的错误： 变量“as.name(var1)”的类型（符号）无效

我怎样才能解决这个问题并让它工作？

Answer 1

您似乎想要使用两个变量的所有组合的模型。这是另一种方法，使用内置的mtcars 数据框进行说明，并使用mpg 作为结果变量。

我们使用combn 获得两个变量的所有组合（不包括结果变量，在本例中为mpg）。 combn 返回一个列表，其中每个列表元素是一个包含一对变量名称的向量。然后我们使用map（来自purrr 包）为每对变量创建模型并将结果存储在一个列表中。

我们使用reformulate 来构造模型公式。 .x 引用变量名称的向量（vars 的每个元素）。例如，如果您运行reformulate(paste(c("cyl", "disp"),collapse="*"), "mpg")，您可以看到reformulate 在做什么。

library(purrr)

# Get all combinations of two variables
vars = combn(names(mtcars)[-grep("mpg", names(mtcars))], 2, simplify=FALSE)

现在我们要对所有变量对运行回归模型并将结果存储在一个列表中：

# No interaction
models = map(vars, ~ glm(reformulate(.x, "mpg"), data=mtcars))

# Interaction only (no main effects)
models = map(vars, ~ glm(reformulate(paste(.x, collapse=":"), "mpg"), data=mtcars))

# Interaction and main effects
models = map(vars, ~ glm(reformulate(paste(.x, collapse="*"), "mpg"), data=mtcars))

用该模型的公式命名每个列表元素：

names(models) = map(models, ~ .x[["terms"]])

要使用paste 而不是reformulate 创建模型公式，您可以这样做（将+ 更改为: 或*，具体取决于您想要包含的交互作用和主要影响的组合）：

models = map(vars, ~ glm(paste("mpg ~", paste(.x, collapse=" + ")), data=mtcars))

要查看这里如何使用paste，您可以运行：

paste("mpg ~", paste(c("cyl", "disp"), collapse=" * "))

当模型同时包含主效应和交互作用时，前两个模型如下所示：

models[1:2]

$`mpg ~ cyl * disp`

Call:  glm(formula = reformulate(paste(.x, collapse = "*"), "mpg"), 
    data = mtcars)

Coefficients:
(Intercept)          cyl         disp     cyl:disp  
   49.03721     -3.40524     -0.14553      0.01585  

Degrees of Freedom: 31 Total (i.e. Null);  28 Residual
Null Deviance:        1126 
Residual Deviance: 198.1  AIC: 159.1

$`mpg ~ cyl * hp`

Call:  glm(formula = reformulate(paste(.x, collapse = "*"), "mpg"), 
    data = mtcars)

Coefficients:
(Intercept)          cyl           hp       cyl:hp  
   50.75121     -4.11914     -0.17068      0.01974  

Degrees of Freedom: 31 Total (i.e. Null);  28 Residual
Null Deviance:        1126 
Residual Deviance: 247.6  AIC: 166.3

要评估模型输出，您可以使用 broom 包中的函数。下面的代码分别返回数据帧，其中包含每个模型的系数和性能统计信息。

library(broom)

model_coefs = map_df(models, tidy, .id="Model")
model_performance = map_df(models, glance, .id="Model")

以下是具有主效应和交互作用的模型的结果：

head(model_coefs, 8)

             Model        term    estimate   std.error statistic      p.value
1 mpg ~ cyl * disp (Intercept) 49.03721186 5.004636297  9.798357 1.506091e-10
2 mpg ~ cyl * disp         cyl -3.40524372 0.840189015 -4.052950 3.645320e-04
3 mpg ~ cyl * disp        disp -0.14552575 0.040002465 -3.637919 1.099280e-03
4 mpg ~ cyl * disp    cyl:disp  0.01585388 0.004947824  3.204212 3.369023e-03
5   mpg ~ cyl * hp (Intercept) 50.75120716 6.511685614  7.793866 1.724224e-08
6   mpg ~ cyl * hp         cyl -4.11913952 0.988229081 -4.168203 2.672495e-04
7   mpg ~ cyl * hp          hp -0.17068010 0.069101555 -2.469989 1.987035e-02
8   mpg ~ cyl * hp      cyl:hp  0.01973741 0.008810871  2.240120 3.320219e-02

Answer 2

您可以按照@akrun 的建议使用fit <- glm(as.formula(paste0("re78 ~ ", var1)), data=newData)。此外，您可能不想调用您的对象glm.fit，因为有一个具有相同功能的函数。

警告：我不明白为什么你有双循环和:。您不想使用单个 covaraite 进行回归吗？我不知道你想达到什么目的。