Испытание anova не выполняется на подгонках, созданных с помощью вставленной формулы - PullRequest
8 голосов
/ 05 октября 2011

Я часто задаю аргумент формулы для функций подгонки модели, таких как lm или lme, вставляя нужные мне части, как в ответе @ DWin на этот вопрос: Понимание lm и среды .

На практике это выглядит так:

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"))

Однако при использовании lme из пакета nlme сравнение двух моделей, построенных таким образом, с использованием anova, не работает, поскольку anova.lme просматривает аргумент сохраненной формулы, чтобы убедиться, что модели были подогнаны под тот же ответ, и сохраненный аргумент формулы просто as.formula(f). Ошибка:

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

Вот что должна делать команда anova (переоборудуя модели, чтобы она работала):

> 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

Есть предложения?

Ответы [ 2 ]

13 голосов
/ 06 октября 2011

Ответ Бена работает, но do.call предоставляет более общее решение, которое он хотел.

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

Он работает, потому что (по умолчанию) аргументы в args = оцениваются перед передачей в lme.

3 голосов
/ 06 октября 2011

Вот взлом, который, кажется, работает:

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
}

но мне это совсем не нравится.Мне бы очень хотелось увидеть более принципиальный ответ на этот вопрос, потому что я постоянно сталкиваюсь с подобной проблемой при попытке расширить пакет bbmle.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...