Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

anova test fails on lme fits created with pasted formula

Tags:

r

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?

like image 958
Aaron left Stack Overflow Avatar asked Oct 05 '11 19:10

Aaron left Stack Overflow


2 Answers

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.

like image 123
Josh O'Brien Avatar answered Nov 13 '22 07:11

Josh O'Brien


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.

like image 4
Ben Bolker Avatar answered Nov 13 '22 09:11

Ben Bolker