在 R 中的重新公式化（）项上使用 anova()答案

【问题标题】：Using anova() on reformulated() terms in R在 R 中的重新公式化（）项上使用 anova()
【发布时间】：2016-08-22 14:37:45
【问题描述】：

背景：

我编写了一个包装函数，允许用户指定数据集中变量的名称，以用作 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 是 Level-2（集群级别）预测变量，x 是 Level-1（单元级别）预测变量，cluster_label 是定义唯一集群的变量名称。例如，如果w = c("w1","w2)、x = "x1" 和我的集群由school 定义，那么（就拟合模型而言）这相当于调用：

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

问题：

函数运行得很好。我的问题是我希望能够使用anova() 来比较两个这样的模型对象，但是当我尝试这样做时，我收到以下错误，由于我使用了reformulate：

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

事实上，当我打电话给summary(mod) 时，我会在其他返回的信息中得到以下行：

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

我注意到anova() 允许通过... 传递额外的参数，那么有没有办法使用这个参数来使函数按预期运行？或者，是否有另一种方法可以对我返回的两个模型对象执行似然比测试而不会出现此错误？如果我需要提供更多详细信息，请告诉我。

可重现的示例：

这些是 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)

【问题讨论】：

您能否提供一个可重现的示例以方便测试？
刚刚添加的例子
啊，我知道有更好的方法来读取该 .csv 文件。感谢编辑！

标签： r model wrapper hierarchical-data nlme

【解决方案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.

【讨论】：