9

バッチで尤度比テストを行うラッパー関数を作成しようとしました。初期モデルを更新するために update() を含めようとしました。ただし、関数内でオブジェクトを探すのではなく、グローバル環境でオブジェクトを探すようです。

fake <- data.frame(subj= rep(1:5, 4), 
                   factor1 = rep(LETTERS[c(1,2,1,2)], each=5), 
                   factor2 = rep(letters[1:2], each=10), 
                   data=sort(rlnorm(20)))

foo <- function(){
                  temp <- fake
                  model1 <- lmer(data~factor1*factor2 + (1 |subj), temp)
                  model1a <- update(model1, ~.-factor1:factor2)
                  model1a}

そして、以下のエラーメッセージが表示されます。

Error in eval(expr, envir, enclos) : object 'factor1' not found

関数内で update() 検索を行う方法はありますか? ありがとうございました！

編集：

私はミスを犯した。「fake」ではなく「temp」をlmerに渡したかったのです。

EDIT2: 提案された便利な解決策の 1 つは、単にデータ オブジェクトを指定することです。update() はこれで問題ありませんが、anova() は、比較しようとしているモデルが異なるデータ オブジェクトに基づいていると考えているようです。

 foo <- function(){
                  temp <- fake
                  model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
                  model1a <- update(model1, ~.-factor1:factor2, data=temp)
                  anova(model1, model1a)
            }
 foo()

エラー メッセージが表示されます。

 Error in anova(model1, model1b) : 
   all models must be fit to the same data object

このエラーは update() を超えていると思います。しかし、これを解決する方法を誰かが知っているのだろうか。update() を使用せずに関数を記述し、代わりにモデルを詳しく説明すると (以下を参照)、上記のエラーは解消されることに注意してください。

 foo <- function(){
                  temp <- fake
                  model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
                  model1a <- lmer(data~factor1 + factor2 + (1 |subj), data=temp)
                  anova(model1, model1a)
            }
 foo()

 Data: temp
 Models:
 model1a: data ~ factor1 + factor2 + (1 | subj)
 model1: data ~ factor1 * factor2 + (1 | subj)
         Df     AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)  
 model1a  5 -4.6909 3.7535  7.3454                           
 model1   6 -8.8005 1.3327 10.4003 6.1097      1    0.01344 *
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

編集 3: 問題は anova() にあるようです。@hadleyの提案も試しました

foo2 <- function(){
  my_update <- function(mod, formula = NULL, data = NULL) {
  call <- getCall(mod)
  if (is.null(call)) {
    stop("Model object does not support updating (no call)", call. = FALSE)
  }
  term <- terms(mod)
  if (is.null(term)) {
    stop("Model object does not support updating (no terms)", call. = FALSE)
  }
  if (!is.null(data)) call$data <- data
  if (!is.null(formula)) call$formula <- update.formula(call$formula, formula)
  env <- attr(term, ".Environment")
  eval(call, env, parent.frame())}

      model1 <- lmer(data~factor1*factor2 + (1 |subj), temp)
      model1a <- my_update(model1, ~.-factor1:factor2)
      anova(model1, model1a)
 }
 foo2()

以下に示すようなエラーメッセージが表示されました。

 Error in as.data.frame.default(data) : 
   cannot coerce class 'structure("mer", package = "lme4")' into a data.frame
4

2 に答える 2

9

私も以前にこの動作に悩まされていたので、独自のバージョンのupdate. 数式の環境ですべてを評価するため、かなり堅牢なはずです。

my_update <- function(mod, formula = NULL, data = NULL) {
  call <- getCall(mod)
  if (is.null(call)) {
    stop("Model object does not support updating (no call)", call. = FALSE)
  }
  term <- terms(mod)
  if (is.null(term)) {
    stop("Model object does not support updating (no terms)", call. = FALSE)
  }

  if (!is.null(data)) call$data <- data
  if (!is.null(formula)) call$formula <- update.formula(call$formula, formula)
  env <- attr(term, ".Environment")

  eval(call, env, parent.frame())
}

library(nlme4)

fake <- data.frame(
  subj = rep(1:5, 4), 
  factor1 = rep(LETTERS[c(1,2,1,2)], each = 5), 
  factor2 = rep(letters[1:2], each = 10), 
  data = sort(rlnorm(20)))

foo <- function() {
  temp <- fake
  model1 <- lmer(data ~ factor1 * factor2 + (1 | subj), fake)
  model1a <- my_update(model1, ~ . - factor1:factor2)
  model1a
}
foo()
于 2012-12-03T20:14:33.550 に答える
4

私は@Hadleyの答えが本当に好きですが（そしておそらく自分でその関数を使用するでしょう）、関数でdata引数を指定することもできますupdate。(ここでは、モデルに渡したいとtemp思います。)

model1a <- update(model1, ~.-factor1:factor2, data = temp)

編集

モデルを と比較しようとしている場合anova、は引数updateの名前を改ざんし、2 つのモデルが 2 つの異なるデータセットを使用して適合したと信じ込ませます。数式のみを更新して新しいモデルを作成すると、これを回避できます。dataanova

foo <- function(){
                  temp <- fake
                  model1 <- lmer(data~factor1*factor2 + (1 |subj), data=temp)
                  newForm <- update.formula(formula(model1), ~.-factor1:factor2)
                  model1a <- lmer(newForm, data=temp)
                  anova(model1, model1a)
            }
于 2012-12-03T20:31:54.947 に答える