사용 R의 재 공식화() 조건에 대한 분산 분석()

Question

배경 : 나는 사용자가 NLME의 lme() 기능에 예측 인자로 사용할 데이터 세트에 변수의 이름을 지정할 수있는 래퍼 함수를 ​​쓴

계층 적 선형 모델에 적합합니다. 관련 부분은 다음과 같습니다

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

곳 w 레벨 2 (클러스터 수준) 예측, x는 레벨 1 (단위 수준) 예측 인자이며, cluster_label 고유 정의하는 변수의 이름입니다 클러스터.

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

문제 :

기능 실행을 w = c("w1","w2)이 x = "x1", 내 클러스터가 (모델을 피팅하는 간다까지로) 다음 school에 의해 정의 된 경우 예를 들어,이 호출하는 것과 같습니다 잘 됐네. ,

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

와 : 내 문제로 인해 reformulate 내 사용 나는이 같은 모델 객체를 비교하는 anova()을 사용할 수 있도록하려면,하지만 나는 다음과 같은 오류가이 작업을 수행 할 때, 내가 수집이다 내가 summary(mod)를 호출 할 때 실제로, 나는 다른 반환 된 정보 중 다음 줄을 얻을 그래서이 인수를 사용하는 방법이있다, 나는 anova() 추가 인수가 ...를 통해 전달 될 수 있습니다 것을 알 수

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

을되도록 기능 원하는대로 동작합니까? 또는이 오류가 발생하지 않고 반환 된 모델 객체 중 두 개에 대해 우도 비율 테스트를 수행하는 다른 방법이 있습니까? 추가 세부 정보를 제공해야하는지 알려주세요.

재현 예 : 이들은 잘 예를 재현 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.