在函数内使用局部协变量更新模型的update()函数

10

我需要在函数内更新一个回归模型。理想情况下,该函数应适用于任何类型的模型(lmglmmultinomclm)。更确切地说,我需要添加一个或多个在函数内定义的协变量。以下为示例:

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)

在评估(expr, envir, enclos)时出错:找不到对象'randData'

这是使用glm的另一个例子

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

有什么想法吗?

2个回答

12

问题在于,在数据框和模型环境中查找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))
}

以上可能是本答案中呈现的最佳解决方案,因为它避免了对内部结构进行操作。它似乎适用于lmglmmultinomclm 。下面的其他解决方案会对内部结构进行修改,因此在模型拟合例程上不太通用。其他解决方案都适用于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),但它是最短的。它们适用于lmglmmultinomclm

请注意,在以下所有示例中,我们实际上并不需要将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) 
}

#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,有时会被人所不齿。它适用于lmglm,但对于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)中的解决方案可以在测试部分运行所有示例。


非常感谢您提供详细和全面的答案!我相信使用环境是正确的方法,但它与glm不兼容(我不明白为什么)。 - Matthias Studer
1
已经进行了修订,使得第一个解决方案适用于 lmglmmultinomclm。它是最通用的,因为它对内部的访问最少。上面的其他解决方案至少适用于 lm - G. Grothendieck
这真的是一个很好的答案。谢谢你解释为什么第一个解决方案是最好的。 - Matthias Studer
已简化(1)中的 MyUpdate - G. Grothendieck
已经改进了第(2)部分的3种解决方案,现在它们可以在测试部分运行4个示例而没有错误。然而,第(3)部分仍然只能运行其中的一部分示例。 - G. Grothendieck

3

一些理论。一个公式对象通常会有一个关联的环境:

frm1 <- y~x # a formula created in the global environment ("in the console")
attr(frm1, ".Environment") # see also unclass(frm1)
## <environment: R_GlobalEnv>

在这里,作用于 frm1 的函数将会知道它们应该在全局环境中寻找 yx(除非另有说明,如 lm() 函数的 data 参数)。另一方面:

f <- function() { y~x }; frm2 <- f() # a formula created in a function
attr(frm2, ".Environment")
## <environment: 0x2f87e48>

这样的公式指向f()中的“局部变量”yx
如果你将在全局环境中创建的公式传递给函数,你通常无法引用在该函数中创建的对象。 解决方案。 通过lm()返回的对象,在其底层包含了相应的公式和环境,但它们是“隐藏”的。但是,可以通过以下代码来解决此问题。
MyUpdate <- function(model){
     assign("randData", data.frame(var1=rnorm(length(model$residuals))),
        envir=attr(model1$terms, ".Environment"))
     model2 <- update(model, ".~.+randData$var1")
     return(model2)
}

感谢您的全面回答!但是,我不想修改原始环境。 - Matthias Studer

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接