使用してRで改質さ（）条件でANOVA（）

Question

背景：私は、ユーザーがNLMEのlme()機能に予測因子として使用するデータセット内の変数の名前を指定することができますラッパー関数を書いた

階層的な線形モデルに適合します。 xは、レベル1（単位レベル）の予測因子であり、cluster_labelが一意の定義変数の名前であり、wは、レベル2（クラスタレベル）である

predictors <- str_c(c(w,x), collapse = " + ") 
mod = lme(fixed = reformulate(termlabels = predictors, response = y), 
      random = reformulate(termlabels = paste0("1|", cluster_label)), 
      data = dat) 

予測子：関連部分はこのようになりますクラスタ。例えば、x = "x1"、w = c("w1","w2)あれば、そして私のクラスタは、これは呼び出すことと同じです（限りモデルが行くフィッティングなど）、その後、schoolによって定義されています。

mod = lme(fixed = y ~ w1 + w2 + x1, random = ~1|school, data = tmp) 

問題：

機能の実行ちょうど良い。私の問題は、私は2つのこのようなモデルオブジェクトを比較するためにanova()を使用できるようにしたいということですが、私は、私は次のエラーを取得するこれを実行しようとすると、私が原因reformulateの私の使用に集まる：

Error in reformulate(termlabels = predictors, response = y) : 
    object 'predictors' not found 

と、私は summary(mod)を呼び出すときに実際に、私は他の返された情報の中で次の行を取得：そう、この引数を使用する方法があり、私は anova()は、追加の引数は ...経由で渡すことができますことに気づく

Fixed: reformulate(termlabels = predictors, response = y) 

をするように機能必要に応じて動作しますか？あるいは、返されたモデルオブジェクトのうちの2つについてこのエラーを出さずに尤度比テストを実行する別の方法がありますか？追加情報を提供する必要がある場合は、私にお知らせください。

再現例：

これらはちょうど良い例を再現HSB12データです。これにはnlmeとstringrのパッケージが必要です。

fitVAM <- function(dat,y,w,x,cluster_label) { 
    library(nlme) 
    library(stringr) 
    predictors <- str_c(c(w,x),collapse=" + ") 
    mod = lme(fixed = reformulate(termlabels = predictors, response = y), 
      random = reformulate(termlabels = paste0("1|", cluster_label)), 
      data = dat) 
    return(mod) 
} 
dat <- read.csv("http://www.ats.ucla.edu/stat/paperexamples/singer/hsb12.csv") 
mod1 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = "female", cluster_label = "school") 
mod2 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = c("female", "ses"), cluster_label = "school") 
anova(mod1,mod2) 

Answer 1

あり、よりストレートな方法かもしれないが、あなたはdo.callを使用することができます。

fitVAM <- function(dat,y,w,x,cluster_label) { 
    library(nlme) 
    library(stringr) 
    predictors <- str_c(c(w,x),collapse=" + ") 
    params <- list(fixed = reformulate(termlabels = predictors, response = y), 
       random = reformulate(termlabels = paste0("1|", cluster_label)), 
       data = dat) 
    mod = do.call(lme, params) 
    mod$call[[3]] <- substitute(dat) #otherwise dat is shown expanded in summary and print output 
    return(mod) 
} 
#dat <- read.csv("http://www.ats.ucla.edu/stat/paperexamples/singer/hsb12.csv") 
mod1 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = "female", cluster_label = "school") 
mod2 <- fitVAM(dat = dat, y = "mathach", w = "meanses", x = c("female", "ses"), cluster_label = "school") 
anova(mod1,mod2) 

#  Model df  AIC  BIC logLik Test L.Ratio p-value 
#mod1  1 5 46908.59 46942.99 -23449.29       
#mod2  2 6 46530.02 46571.29 -23259.01 1 vs 2 380.5734 <.0001 
#Warning message: 
#In anova.lme(mod1, mod2) : 
# fitted objects with different fixed effects. REML comparisons are not meaningful.