自动公式构建答案

【问题标题】：Automated formula construction自动公式构建
【发布时间】：2013-08-30 17:25:21
【问题描述】：

在最近的一项家庭作业中，我们被要求运行 27 个线性模型，每次都添加一个附加变量（目标是绘制 R² 的变化与调整后的 R²）。我发现很难通过算法创建这样的公式。我最终使用的代码如下所示（请注意，数据框中的第一列是因变量，其余都是预期的自变量。

 make.formula <- function(howfar) {
  formula <- c()
  for (i in 1:howfar) {
    if (i == 1) {
      formula <- paste(formula, names(d)[i], '~')}
    else if (i == howfar) {
      formula <- paste(formula, names(d)[i], '')
    }
    else {
      formula <- paste(formula, names(d)[i], '+')}
  }
  return(formula)
}

formulas <- lapply(seq(2, length(d)), make.formula)
formulas <- lapply(formulas, as.formula)
fits <- lapply(formulas, lm, data = d)

这可行，但似乎远非理想，我的印象是我在 R 中使用 for 循环所做的任何事情都可能不是最好的方式。有没有更简单的方法来为给定的数据框通过算法构造公式？

【问题讨论】：

标签： r modeling

【解决方案1】：

reformulate()，一个用于从字符向量创建公式的漂亮函数，可能会派上用场。以下是它的作用示例：

reformulate(response="Y", termlabels=c("X1", "X2", "X3"))
# Y ~ X1 + X2 + X3

以下是您在实践中可能会如何使用它。（请注意，我在这里在 lm() 调用内部创建公式。因为 formula 对象带有有关它们创建环境的信息，所以我会有点犹豫在外部创建它们 的 lm() 调用中，您实际上想在其中使用它们。）：

evars <- names(mtcars)[2:5]
ii <- lapply(1:4, seq_len)

lapply(ii, 
       function(X) {
          coef(lm(reformulate(response="mpg", termlabels=evars[X]), data=mtcars))
})
# [[1]]
# (Intercept)         cyl 
#    37.88458    -2.87579 
# 
# [[2]]
# (Intercept)         cyl        disp 
# 34.66099474 -1.58727681 -0.02058363 
# 
# [[3]]
# (Intercept)         cyl        disp          hp 
# 34.18491917 -1.22741994 -0.01883809 -0.01467933 
# 
# [[4]]
# (Intercept)         cyl        disp          hp        drat 
# 23.98524441 -0.81402201 -0.01389625 -0.02317068  2.15404553

【讨论】：

(+1) 使用lapply(1:4, seq_len) 真的很酷。我无法弄清楚如何在没有循环的情况下生成这些字符串。而且我学会了reformulate 函数。谢谢！
乔希的精彩回答。 lapply(1:4, seq_len) 本身就是一个方便的提示，奇怪的是我还没有遇到reformulate 的任何文档。欣赏有关在 lm 调用中指定公式的最佳实践提示。

【解决方案2】：

Map 可以用来解决这个问题：

mydata<-mtcars

dep<-as.list(rep("mpg~",(dim(mydata)[2]-1))) # ldependent variables with ~ 
indep1<- as.list( names(mydata)[-1])
indeno<-as.list(1:(dim(mydata)[2]-1))

myreg<-Map(function(x,y) (lm(as.formula(paste(x,paste(unlist(indep[1:y]),collapse="+"))),data=mtcars))$coefficient,dep,indeno)
   > myreg
[[1]]
(Intercept)         cyl 
   37.88458    -2.87579 

[[2]]
(Intercept)         cyl        disp 
34.66099474 -1.58727681 -0.02058363 

[[3]]
(Intercept)         cyl        disp          hp 
34.18491917 -1.22741994 -0.01883809 -0.01467933 

[[4]]
(Intercept)         cyl        disp          hp        drat 
23.98524441 -0.81402201 -0.01389625 -0.02317068  2.15404553 

[[5]]
(Intercept)         cyl        disp          hp        drat          wt 
36.00835689 -1.10748650  0.01235733 -0.02401743  0.95220742 -3.67328708 

[[6]]
(Intercept)         cyl        disp          hp        drat          wt        qsec 
26.30735899 -0.81856023  0.01320490 -0.01792993  1.32040573 -4.19083238  0.40146117 

[[7]]
(Intercept)         cyl        disp          hp        drat          wt        qsec          vs 
25.88354175 -0.85665309  0.01314097 -0.01733070  1.31265550 -4.22434351  0.44873351 -0.27816899 

[[8]]
(Intercept)         cyl        disp          hp        drat          wt        qsec          vs          am 
15.57313068 -0.27859352  0.01471012 -0.02144242  0.81505862 -3.94373934  0.80975689  0.36835866  2.79374984 

[[9]]
(Intercept)         cyl        disp          hp        drat          wt        qsec          vs          am 
12.83083549 -0.16881263  0.01623358 -0.02424055  0.70590083 -4.03214213  0.86828517  0.36470431  2.55092849 
       gear 
 0.50293618 

[[10]]
(Intercept)         cyl        disp          hp        drat          wt        qsec          vs          am 
12.30337416 -0.11144048  0.01333524 -0.02148212  0.78711097 -3.71530393  0.82104075  0.31776281  2.52022689 
       gear        carb 
 0.65541302 -0.19941925

【讨论】：