无需字符串操作，以编程方式构建公式

Question

无需字符串操作，以编程方式构建公式

4

为了举例说明，在R中考虑一个基本的回归模型：

form1 <- Petal.Length ~ Sepal.Length + Sepal.Width
fit1 <- lm(form1, iris)

（对任何在此发帖的植物学家表示歉意。）

为了添加二次项和交互项，我知道有三种方法：

1）老派的方式

逐个输入项：

form2 <- . ~ Sepal.Length*Sepal.Width + I(Sepal.Length^2) + I(Sepal.Width^2)
fit2 <- update(fit1, form2)

这种方法不适用于复杂的公式，也不能用来编程。

2) 糟糕的方法

字符串操作：

vars <- attr(terms(form1), "term.labels")
squared_terms <- sprintf("I(%s^2)", vars)
inter_terms <- combn(vars, 2, paste, collapse = "*")
form2 <- reformulate(c(inter_terms, squared_terms), ".")

这个可以扩展，但它并不真正可编程，因为函数本身需要硬编码。

3) 后门

直接操作数据。

library(lazyeval)
library(dplyr)

square <- function (v) interp(~ I(v1^2), v1 = as.name(v))
inter <- function(v) interp(~ v1*v2, v1 = as.name(v[1]), v2 = as.name(v[2]))

vars <- attr(terms(form1), "term.labels")
squared_terms <- lapply(vars, square) %>%
  set_names(paste0(vars, " ^2"))
inter_terms <- combn(vars, 2, inter, simplify = FALSE) %>%
  set_names(combn(vars, 2, paste, collapse = " x "))

fit2 <- model.frame(fit1) %>%
  mutate_(.dots = squared_terms) %>%
  mutate_(.dots = inter_terms) %>%
  lm(Petal.Length ~ ., data = .)

这是相当可扩展的，可以编程到变量命名。但它也有点疯狂，因为它打破了首先使用公式的目的。

我希望我能做到的

我希望我能做类似于这样的事情：

library(lazyeval)
library(dplyr)

square <- function (v) interp(~ I(v1^2), v1 = as.name(v))
inter <- function(v) interp(~ v1*v2, v1 = as.name(v[1]), v2 = as.name(v[2]))

squared_terms <- apply.formula(form1, squared_terms)
inter_terms <- combn.formula(form1, 2, inter)

fit2 <- form1 %>%
  append.formula(squared_terms) %>%
  append.formula(inter_terms) %>%
  update(fit1, .)

除了滥用dplyr之外，这里有两个杀手级功能：

1.能够以编程方式从基本的R对象生成任意公式项

2.能够将项添加到公式中，这些项将像手动输入的项一样工作

使用方法3可以获得功能1，使用方法2可以获得功能2。是否有一种方法4 -- "中庸之道" -- 可以同时获取这两种功能？

- shadowtalker

创建设计矩阵 --> lm.fit - rawr

1

您可以轻松地使用 lm(Petal.Length ~ (Sepal.Length + Sepal.Width)^2, iris) 进行所有交互。我猜方括号里的平方部分可能会有些棘手。 - MrFlick

@rawr 这就是我担心的事情。 - shadowtalker

为什么？计算并将其放入矩阵比进行字符串和公式操作更容易，适用于 lm。 - rawr

@rawr 整个重点是将所有内容封装在几个简单的函数中。设计矩阵方法涉及大量样板代码。这个想法是为了能够更加灵活地拟合模型。我喜欢将编程和建模分开。 - shadowtalker

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- MrFlick · Accepted Answer

好的，这里涉及到很多方面，但是以下是一些帮助函数，可以执行非常具体的任务。

extract_rhs_symbols <- function(x) {
    as.list(attr(delete.response(terms(x)), "variables"))[-1]
}
symbols_to_formula <- function(x) {
    as.call(list(quote(`~`), x))    
}
sum_symbols <- function(...) {
    Reduce(function(a,b) bquote(.(a)+.(b)), do.call(`c`, list(...), quote=T))
}
square_terms <- function(x) {
    symbols_to_formula(sum_symbols(sapply(extract_rhs_symbols(x), function(x) bquote(I(.(x)^2)))))
}
interact_rhs<-function(x) {
    x[[length(as.list(x))]] <- bquote((.(x[[length(as.list(x))]]))^2)
    x
}
add_rhs_dot <- function(x) {
   x[[length(as.list(x))]] <- sum_symbols(quote(.), x[[length(as.list(x))]])    
   x
}
add_terms<-function(f, x) {
    update(f, add_rhs_dot(x))
}

所有这些基本上都是通过调用公式来操作的。

所以，如果您有一个像下面这样的公式：

my.formula <- Petal.Length ~ Sepal.Length + Sepal.Width + Other

您可以使用以下方法添加平方项：

add_terms(my.formula, square_terms(my.formula))

你可以通过右键单击获取所有的右侧交互。

interact_rhs(my.formula)

或者使用以下方式同时完成

add_terms(interact_rhs(my.formula), square_terms(my.formula))

这提供了

Petal.Length ~ Sepal.Length + Sepal.Width + Other + I(Sepal.Length^2) + 
    I(Sepal.Width^2) + I(Other^2) + Sepal.Length:Sepal.Width + 
    Sepal.Length:Other + Sepal.Width:Other

我还没有进行全面测试，因此可能会存在一些情况导致它失效，但在大多数简单情况下应该可以工作。