【问题标题】:update() a model inside a function with local covariateupdate() 具有局部协变量的函数内的模型
【发布时间】:2014-06-20 03:12:03
【问题描述】:

我需要从函数内部更新回归模型。理想情况下,该函数应该适用于任何类型的模型(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)

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

这是 glm 的另一个例子

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

有什么想法吗?

【问题讨论】:

    标签: r glm lm


    【解决方案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))
    }
    

    以上可能是此答案中提出的最佳解决方案,因为它避免了内部结构的混乱。它似乎适用于lmglmmultinomclm。下面的其他解决方案确实与内部结构混在一起,因此在模型拟合例程中不太通用。其他人都使用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),但它是最短的。他们与lmglmmultinomclm 合作。

    请注意,我们实际上不需要将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,这有时是不受欢迎的。它适用于lmglmclm,除了输出不显示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(我不明白为什么)。
    • 已修改,以便第一个解决方案适用于 lmglmmultinomclm。它是最通用的,因为它对内部的访问最少。上述其他解决方案至少在 lm 上有效。
    • 这真是一个很好的答案。感谢您解释为什么第一个解决方案是最好的解决方案。
    • 在 (1) 中简化了MyUpdate
    • 改进了 (2) 中的 3 个解决方案,现在他们运行测试部分中的 4 个示例而不会出错。 (3) 仍然只运行其中的一个子集。
    【解决方案2】:

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

    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 arg)。另一方面:

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

    这样的公式指向yxf() 中的“局部变量”。

    如果您将在全局环境中创建的公式传递给函数,则大多数情况下您将无法引用在该函数中创建的对象。

    解决方案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)
    }
    

    【讨论】:

    • 感谢您的全面回答!但是,我不想修改原来的环境。
    猜你喜欢
    • 1970-01-01
    • 2022-08-02
    • 2016-12-30
    • 1970-01-01
    • 2020-04-04
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多