I tried to write a wrapper function to do likelihood ratio tests in batches. I tried to include update() to update the initial model. However, it seems that instead of looking for objects inside the function, it searches for objects in the global environment.
我试着编写一个包装函数来批量进行似然比检验。我试图包含update()来更新初始模型。但是,它似乎不是在函数内部查找对象,而是在全局环境中搜索对象。
fake <- data.frame(subj= rep(1:5, 4),
factor1 = rep(LETTERS[c(1,2,1,2)], each=5),
factor2 = rep(letters[1:2], each=10),
data=sort(rlnorm(20)))
foo <- function(){
temp <- fake
model1 <- lmer(data~factor1*factor2 + (1 |subj), temp)
model1a <- update(model1, ~.-factor1:factor2)
model1a}
And it gives an error message below:
它在下面给出了一条错误消息:
Error in eval(expr, envir, enclos) : object 'factor1' not found
Is there anyway to make update() search within the function? Thank you!
无论如何在函数内进行update()搜索?谢谢!
EDIT:
编辑:
I made a mistake. I wanted to pass "temp" to lmer, not "fake".
我犯了一个错误。我想把“temp”传给lmer,而不是“假”。
EDIT2: One convenient solution suggested is to simply specify the data object. Although update() now has no problem with this, anova() seems to think that the models I am trying to compare are based on different data objects
EDIT2:建议的一个方便的解决方案是简单地指定数据对象。虽然update()现在没有问题,但anova()似乎认为我试图比较的模型是基于不同的数据对象
foo <- function(){
temp <- fake
model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
model1a <- update(model1, ~.-factor1:factor2, data=temp)
anova(model1, model1a)
}
foo()
I get an error message:
我收到一条错误消息:
Error in anova(model1, model1b) :
all models must be fit to the same data object
I suppose this error goes beyond update(). But I wonder if anyone knows how this can be resolved. Note that if I write the function without using update() and instead spell out the models (see below), the error above goes away:
我想这个错误超出了update()。但我想知道是否有人知道如何解决这个问题。请注意,如果我在不使用update()的情况下编写函数,而是拼出模型(见下文),则上面的错误就会消失:
foo <- function(){
temp <- fake
model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
model1a <- lmer(data~factor1 + factor2 + (1 |subj), data=temp)
anova(model1, model1a)
}
foo()
Data: temp
Models:
model1a: data ~ factor1 + factor2 + (1 | subj)
model1: data ~ factor1 * factor2 + (1 | subj)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model1a 5 -4.6909 3.7535 7.3454
model1 6 -8.8005 1.3327 10.4003 6.1097 1 0.01344 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
EDIT 3: It seems that the issue is with anova(). I also tried the suggestion by @hadley
编辑3:似乎问题出在anova()上。我也试过@hadley的建议
foo2 <- function(){
my_update <- function(mod, formula = NULL, data = NULL) {
call <- getCall(mod)
if (is.null(call)) {
stop("Model object does not support updating (no call)", call. = FALSE)
}
term <- terms(mod)
if (is.null(term)) {
stop("Model object does not support updating (no terms)", call. = FALSE)
}
if (!is.null(data)) call$data <- data
if (!is.null(formula)) call$formula <- update.formula(call$formula, formula)
env <- attr(term, ".Environment")
eval(call, env, parent.frame())}
model1 <- lmer(data~factor1*factor2 + (1 |subj), temp)
model1a <- my_update(model1, ~.-factor1:factor2)
anova(model1, model1a)
}
foo2()
I got an error message as shown below:
我收到如下所示的错误消息:
Error in as.data.frame.default(data) :
cannot coerce class 'structure("mer", package = "lme4")' into a data.frame
2 个解决方案
#1
8
I've been bitten by this behaviour before too, so I wrote my own version of update
. It evaluates everything in the environment of the formula, so it should be fairly robust.
我之前也被这种行为所困扰,所以我编写了自己的更新版本。它评估公式环境中的所有内容,因此它应该相当健壮。
my_update <- function(mod, formula = NULL, data = NULL) {
call <- getCall(mod)
if (is.null(call)) {
stop("Model object does not support updating (no call)", call. = FALSE)
}
term <- terms(mod)
if (is.null(term)) {
stop("Model object does not support updating (no terms)", call. = FALSE)
}
if (!is.null(data)) call$data <- data
if (!is.null(formula)) call$formula <- update.formula(call$formula, formula)
env <- attr(term, ".Environment")
eval(call, env, parent.frame())
}
library(nlme4)
fake <- data.frame(
subj = rep(1:5, 4),
factor1 = rep(LETTERS[c(1,2,1,2)], each = 5),
factor2 = rep(letters[1:2], each = 10),
data = sort(rlnorm(20)))
foo <- function() {
temp <- fake
model1 <- lmer(data ~ factor1 * factor2 + (1 | subj), fake)
model1a <- my_update(model1, ~ . - factor1:factor2)
model1a
}
foo()
#2
4
Although I really like @Hadley's answer (and will likely use that function myself), you can also specify a data
argument in the update
function. (Here, I assumed you wanted to pass temp
to your models.)
虽然我非常喜欢@Hadley的答案(并且可能会自己使用该函数),但您也可以在更新函数中指定数据参数。 (在这里,我以为你想把温度传递给你的模特。)
model1a <- update(model1, ~.-factor1:factor2, data = temp)
EDIT
编辑
If you're looking to compare models with anova
, update
will mung up the name of the data
argument and "trick" anova
into believing that the two models were fit using two different datasets. Updating only the formula and creating a new model will avoid this:
如果您希望将模型与anova进行比较,则更新将清除数据参数的名称并“欺骗”anova,使其相信这两个模型适合使用两个不同的数据集。仅更新公式并创建新模型将避免这种情况:
foo <- function(){
temp <- fake
model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
newForm <- update.formula(formula(model1), ~.-factor1:factor2)
model1a <- lmer(newForm, data=temp)
anova(model1, model1a)
}
#1
8
I've been bitten by this behaviour before too, so I wrote my own version of update
. It evaluates everything in the environment of the formula, so it should be fairly robust.
我之前也被这种行为所困扰,所以我编写了自己的更新版本。它评估公式环境中的所有内容,因此它应该相当健壮。
my_update <- function(mod, formula = NULL, data = NULL) {
call <- getCall(mod)
if (is.null(call)) {
stop("Model object does not support updating (no call)", call. = FALSE)
}
term <- terms(mod)
if (is.null(term)) {
stop("Model object does not support updating (no terms)", call. = FALSE)
}
if (!is.null(data)) call$data <- data
if (!is.null(formula)) call$formula <- update.formula(call$formula, formula)
env <- attr(term, ".Environment")
eval(call, env, parent.frame())
}
library(nlme4)
fake <- data.frame(
subj = rep(1:5, 4),
factor1 = rep(LETTERS[c(1,2,1,2)], each = 5),
factor2 = rep(letters[1:2], each = 10),
data = sort(rlnorm(20)))
foo <- function() {
temp <- fake
model1 <- lmer(data ~ factor1 * factor2 + (1 | subj), fake)
model1a <- my_update(model1, ~ . - factor1:factor2)
model1a
}
foo()
#2
4
Although I really like @Hadley's answer (and will likely use that function myself), you can also specify a data
argument in the update
function. (Here, I assumed you wanted to pass temp
to your models.)
虽然我非常喜欢@Hadley的答案(并且可能会自己使用该函数),但您也可以在更新函数中指定数据参数。 (在这里,我以为你想把温度传递给你的模特。)
model1a <- update(model1, ~.-factor1:factor2, data = temp)
EDIT
编辑
If you're looking to compare models with anova
, update
will mung up the name of the data
argument and "trick" anova
into believing that the two models were fit using two different datasets. Updating only the formula and creating a new model will avoid this:
如果您希望将模型与anova进行比较,则更新将清除数据参数的名称并“欺骗”anova,使其相信这两个模型适合使用两个不同的数据集。仅更新公式并创建新模型将避免这种情况:
foo <- function(){
temp <- fake
model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
newForm <- update.formula(formula(model1), ~.-factor1:factor2)
model1a <- lmer(newForm, data=temp)
anova(model1, model1a)
}