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但可能不适合其他人.

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)以下内容修改原始模型的环境以指向包含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)
}

有关更多信息,请阅读本文档中的代理部分:http://r-proto.googlecode.com/files/prototype_approaches.pdf

2b)这可以在没有原型的情况下完成,但代价是将##行扩展到包含一些额外丑陋环境操作的三行.这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)中获得解决方案以运行测试部分中的所有示例.