anova test fails on lme fits created with pasted formula

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.

anova test fails on lme fits created with pasted formula

Related

Recent Posts