问题是 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,但可能不适用于其他人。
test 这是一个测试,如果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) 环境混乱
此外,这里还有三个涉及弄乱环境的解决方案。第一个是最干净的,其次是第二个,然后是第三个。第三个是最不可接受的,因为它实际上将var1 写入模型的环境(危险地覆盖了那里的任何var1),但它是最短的。他们与lm、glm、multinom 和clm 合作。
请注意,我们实际上不需要将var1 放入数据框中,也没有必要将更新公式放在引号中,我们在下面的所有示例中都进行了更改。 return 语句也可以删除,我们也这样做了。
2a) 下面修改原始模型的环境,指向一个新的代理proto对象,该对象包含var1,其父对象是原始模型环境。这里proto(p, var1 = rnorm(n))是代理proto对象(proto对象是一个语义不同的环境),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) 最短但最老套的方法是将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/substitute 这个解决方案确实使用了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) 中获得了解决方案,以运行测试部分中的所有示例。