r - update() 具有局部协变量的函数内的模型

Question

我需要从函数内部更新回归模型。理想情况下，该函数应该适用于任何类型的模型(lm、glm、multinom、clm)。更准确地说，我需要添加一个或多个在函数内部定义的协变量。这是一个例子。

MyUpdate <- function(model){
     randData <- data.frame(var1=rnorm(length(model$residuals)))
     model2 <- update(model, ".~.+randData$var1")
     return(model2)
}

这是一个使用示例

data(iris)
model1 <- lm(Sepal.Length~Species, data=iris)
model2 <- MyUpdate(model1)

eval 中的错误(expr、envir、enclos):找不到对象“randData”

这是 glm 的另一个例子

model1 <- glm(Sepal.Length>5~Species, data=iris, family=binomial)
model2 <- MyUpdate(model1)

任何的想法？

Answer 1

问题是 var1在数据框和模型的环境中查找，但不在 MyUpdate 中的环境中查找.
1) 为避免此问题，不仅要使用修改后的公式更新模型，还要使用包含 var1 的修改后数据框。 :

MyUpdate <- function(model) {
     mf <- model.frame(model)
     n <- nrow(mf)
     var1 <- rnorm(n)
     update(model, formula = . ~ . + var1, data = data.frame(mf, var1))
}

以上可能是此答案中提出的最佳解决方案，因为它避免了内部结构的困惑。它似乎适用于 lm , glm , multinom和 clm .下面的其他解决方案确实与内部结构混在一起，因此在模型拟合例程中不太通用。其他人都使用 lm但可能不适用于其他人。
测试这是一个测试，如果 MyUpdate 在问题中提到的每个模型拟合函数上运行而没有错误如上所述，(2)中的解决方案都运行测试没有错误。解决方案 (3) 至少适用于 lm .

model.lm <- lm(Sepal.Length~Species, data=iris)
MyUpdate(model.lm)

model.glm <- glm(Sepal.Length~Species, data=iris)
MyUpdate(model.glm)

library(nnet)
example(multinom)
MyUpdate(bwt.mu)

library(ordinal)
model.clm <- clm(rating ~ temp * contact, data = wine)
MyUpdate(model.clm)

其余的解决方案执行更直接的内部访问，使其对更改模型函数的鲁棒性降低。
2) 环境困惑
此外，这里还有三个涉及环境困惑的解决方案。第一个是最干净的，其次是第二个，然后是第三个。第三个是最 Not Acceptable ，因为它实际上写 var1进入模型的环境(危险地覆盖任何 var1)，但它是最短的。他们与 lm 合作, glm multinom和 clm .
请注意，我们实际上并不需要输入 var1到数据框中，也没有必要将更新公式放在引号中，我们在下面的所有示例中都进行了更改。还有 return声明可以删除，我们也这样做了。
2a) 下面修改原模型的环境，指向一个新的代理proto对象，其中包含 var1其父级是原始模型环境。这里 proto(p, var1 = rnorm(n))是代理原型(prototype)对象(原型(prototype)对象是具有不同语义的环境)和 p是代理的父级。

library(proto)

MyUpdate <- function(model){

     mf <- model.frame(model)
     n <- nrow(mf)
     var1 <- rnorm(n)
     p <- environment(formula(model))

     if (is.null(model$formula)) {
           attr(model$terms, ".Environment") <- proto(p, var1 = var1)
     } else environment(model$formula) <- proto(p, var1 = var1)

     update(model, . ~ . + var1) 
}

#note: the period is shorthand for 
keep everything on either the left or right hand side 
of the formula (i.e., the ~) and 
the + or - sign are used to add or remove model terms

如需更多信息，请阅读本文档中的代理部分: http://r-proto.googlecode.com/files/prototype_approaches.pdf
2b) 这可以在没有 proto 的情况下交替完成，但代价是将 ## 行扩展到包含一些额外的丑陋环境操作的三行。这里 e是代理环境。

MyUpdate <- function(model){
     mf <- model.frame(model)
     n <- nrow(mf)
     var1 <- rnorm(n)
     p <- environment(formula(model))

     e <- new.env(parent = p)
     e$var1 <- var1

     if (is.null(model$formula)) attr(model$terms, ".Environment") <- e
     else environment(model$formula) <- e

     update(model, . ~ . + var1)
}

2c) 最短但最 hackish 是坚持 var1成原 model环境:

MyUpdate <- function(model){
     mf <- model.frame(model)
     n <- nrow(mf)
     var1 <- rnorm(n)       

     if (is.null(model$formula)) attr(model$terms, ".Environment")$var1 <- var1
     else environment(model$formula)$var1 <- var1

     update(model, . ~ . + var1)
}

3) 评估/替代此解决方案确实使用 eval这有时是不受欢迎的。它适用于 lm和 glm和 clm它可以工作，只是输出不显示 var1而是计算它的表达式。

MyUpdate <- function(model) {
     m <- eval.parent(substitute(update(model, . ~ . + rnorm(nrow(model.frame(model))))))
     m$call$formula <- update(formula(model), . ~ . + var1)
     names(m$coefficients)[length(m$coefficient)] <- "var1"
     m
}

修订 添加了其他解决方案，简化了 (1)，在 (2) 中获得了解决方案，以运行测试部分中的所有示例。