开发者

anova test fails on lme fits created with pasted formula

开发者 https://www.devze.com 2023-04-11 04:59 出处:网络
I often specify the formula argument to model fitting functions like lm or lme by pasting together the parts I need, as in @DWin\'s answer to this question: Understanding lm and environment.

I often specify the formula argument to model fitting functions like lm or lme by pasting together the parts I need, as in @DWin's answer to this question: Understanding lm and environment.

In practice this looks like this:

library(nlme)
set.seed(5)
ns <- 5; ni <- 5; N <- ns*ni
d <- data.frame(y=rnorm(N),
                x1=rnorm(N),
    开发者_如何转开发            x2=factor(rep(1:ni, each=ns)),
                id=factor(rep(1:ns, ni)))

getm <- function(xs) {
  f <- paste("y ~", paste(xs, collapse="+"))
  lme(as.formula(f), random=~1|id, data=d, method="ML")
}
m1 <- getm("x1")
m2 <- getm(c("x1", "x2"))

However, with lme from the nlme package, comparing two models constructed in the way using anova doesn't work, because anova.lme looks at the saved formula argument to ensure that the models were fit on the same response, and the saved formula argument is simply as.formula(f). The error is:

> anova(m1, m2)
Error in inherits(object, "formula") : object 'f' not found

Here's what the anova command should do (refitting the models so that it works):

> m1 <- lme(y~x1, random=~1|id, data=d, method="ML")
> m2 <- lme(y~x1+x2, random=~1|id, data=d, method="ML")
> anova(m1, m2)
   Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m1     1  4 76.83117 81.70667 -34.41558                        
m2     2  8 72.69195 82.44295 -28.34597 1 vs 2 12.13922  0.0163

Any suggestions?


Ben's answer works, but do.call provides the more general solution he wished for.

getm <- function(xs) {
    f <- as.formula(paste("y ~", paste(xs, collapse="+")))
    do.call("lme", args = list(f, random=~1|id, data=d, method="ML"))
}

It works because (by default) the arguments in args = are evaluated before being passed to lme.


Here's a hack that seems to work:

getm <- function(xs) {
  f <- paste("y ~", paste(xs, collapse="+"))
  m <- lme(as.formula(f), random=~1|id, data=d, method="ML")
  m$call$fixed <- eval(m$call$fixed)
  m
}

but I don't like it at all. I would very much like to see a more principled answer to this question, because I run into this sort of problem all the time when trying to extend the bbmle package.

0

精彩评论

暂无评论...
验证码 换一张
取 消

关注公众号