3

拥有一个lm对象，我需要根据其表示为字符向量的变量创建一个函数。我尝试使用eval和的组合expr来创建一个f函数，该函数将进一步用于后者obj并对其进行nlm优化。

library(tidyverse)
df <- drop_na(airquality)
model <- lm(Ozone~. - Temp, data = df, x=TRUE, y=TRUE)
base_vars <- all.vars(formula(model)[-2])
k <- length(base_vars)

f <- function(base_df, x, y, parms) {
  with(base_df, parms[1] + 
         eval(expr(paste(paste(paste0('parms[', 2:(k+1), ']'), base_vars, sep = '*'), collapse = '+'))) + 
         log(parms[k+2] * (x - parms[k+3] ^ 2)))
}
obj <- function(parms, y, x) mean((residuals(model) - f(df, x, y, parms))^2) 
fit <- with(data, nlm(obj, c(0, 0, 0, 0, 0, 0, 0), y = e, x = x))

但是调用f(model$x, df$Temp, model$y, c(0, 0, 0, 0, 0, 0, 0))会导致以下错误：

Error in eval(substitute(expr), data, enclos = parent.frame()) : 
  numeric 'envir' arg not of length one 
4.
eval(substitute(expr), data, enclos = parent.frame()) 
3.
with.default(base_df, parms[1] + eval(expr(paste(paste(paste0("parms[", 
    2:(k + 1), "]"), base_vars, sep = "*"), collapse = "+"))) + 
    log(parms[k + 2] * (x - parms[k + 3]^2))) 
2.
with(base_df, parms[1] + eval(expr(paste(paste(paste0("parms[", 
    2:(k + 1), "]"), base_vars, sep = "*"), collapse = "+"))) + 
    log(parms[k + 2] * (x - parms[k + 3]^2))) 
1.
f(model$x, df$Temp, model$y, c(0, 0, 0, 0, 0, 0, 0))

eval我相信环境与函数所暗示的环境之间可能存在冲突with，但无法弄清楚原因。任何想法如何f为变量模型创建自定义函数？

的预期输出为f(model$x, df$Temp, model$y, c(0, 0, 0, 0, 0, 0, 0))：

with(base_df, parms[1]+parms[2]*Solar.R+parms[3]*Wind+parms[4]*Temp+parms[5]*Month+
              parms[6]*Day+log(parms[7] * (Temp - parms[8] ^ 2)))

但对于不同的模型，它可能是这样的：

with(base_df, 
     parms[1]+parms[2]*var1+parms[3]*var2+log(parms[4]*(var3-parms[5]^2)))

所以每次调用的变量和参数的数量都是不同的。

4

2 回答 2

5

R 支持在该语言上进行计算，但它不应该是您的首选。如果你这样做，它不应该涉及代码的文本处理。您在这里没有需要计算语言的情况。我不知道你认为你的尝试会如何工作，但我不知道这个expr功能，我拒绝安装包 tidyverse 及其巨大的依赖树。

此外，您通常应该避免with在交互式使用之外。但with问题不在这里。

这是我将如何做到这一点：

df <- airquality[complete.cases(airquality),]
model <- lm(Ozone~. - Temp, data = df)

f <- function(base_df, x, parms) {

  m <- model.matrix(model, data = base_df)
  k <- ncol(m)
  stopifnot(length(parms) == (k + 2L))
  #I use exp(parms[k+1]) to ensure a positive value within the log
  m %*% parms[seq_len(k)] + log(exp(parms[k + 1L]) * (x - parms[k + 2L] ^ 2))

}

obj <- function(parms, y, x, base_df) mean((residuals(model) - f(base_df, x, parms))^2) 

#some x:
x <- rpois(nrow(df), 10)

fit <- nlm(obj, c(0, 0, 0, 0, 0, 0, 0), x = x, base_df = df)
#works

您似乎没有使用y，因此我将其从代码中删除。

请注意我如何为线性部分创建设计矩阵（使用model.matrix）并将矩阵乘法与参数一起使用。您还需要确保log不返回Inf// -Inf。NaN

于 2020-02-21T07:23:04.910 回答
1

我认为@Roland 给出了一个很好的答案，涵盖了您的实际问题。我正在根据问题标题来隔离我认为您具体提出的问题，而没有评论这是否是一个好主意。它可能不在这个用例中。

但您所寻找的很可能eval_tidy()来自rlang. 我把::函数符号留在了这里，很明显这里使用的是什么包。

请注意，我修复了一些似乎是代码错误的事情。parms由于日志，我也使用全1而不是零来测试。

library(rlang)
library(tidyr)

# dropped y since it was an unused argument
f <- function(base_df, x, parms) {
  # set an expression to evaluate using parse_expr()
  .f <- rlang::parse_expr(paste(paste(paste0('parms[', 2:(k+1), ']'),
                                      base_vars, sep = '*'), collapse = '+'))

  # use eval_tidy() with the data mask  
  y_part1 <- rlang::eval_tidy(.f, data = base_df)
  y_part2 <- log(parms[k + 2] * (x - parms[k + 3] ^ 2))

  parms[1] + y_part1 + y_part2
}

# using your code
df <- tidyr::drop_na(airquality)
model <- lm(Ozone~. - Temp, data = df, x=TRUE, y=TRUE)
base_vars <- all.vars(formula(model)[-2])
k <- length(base_vars)

# changed to all ones, I think this is what you wanted for length
parms <- rep(1, k + 3)

method_1 <- f(df, df$Temp, parms)

method_2 <- with(df, parms[1]+parms[2]*Solar.R+parms[3]*Wind+parms[4]*Temp+parms[5]*Month+
                   parms[6]*Day+log(parms[7] * (Temp - parms[8] ^ 2)))


all.equal(method_1, method_2)
# [1] TRUE
于 2020-03-07T19:24:05.397 回答