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.
精彩评论