问题在于,在数据框和模型环境中查找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)操作环境
此外还有三个涉及操作环境的解决方案。第一个最干净,其次是第二个,然后是第三个。第三个最不可取,因为它实际上将var1
写入模型的环境中(危险地覆盖任何var1
),但它是最短的。它们适用于lm
、glm
、multinom
和clm
。
请注意,在以下所有示例中,我们实际上并不需要将var1
放入数据框中,也不必将更新公式用引号括起来。同时,return
语句可以被删除,我们也已经这样做了。
2a) 下面的代码修改了原始模型的环境,将其指向一个包含var1
的新代理proto对象,其父对象是原始模型环境。在这里,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)
}
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
,而是显示计算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)中的解决方案可以在测试部分运行所有示例。
lm
、glm
、multinom
和clm
。它是最通用的,因为它对内部的访问最少。上面的其他解决方案至少适用于lm
。 - G. GrothendieckMyUpdate
。 - G. Grothendieck