Using anova() on reformulated() terms in R

Question

Background:

I wrote a wrapper function that allows the user to specify the names of variables in a data set to use as predictors in NLME's lme() function to fit a hierarchical linear model. The relevant part looks like this:

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

where w is Level-2 (cluster-level) predictors, x is Level-1 (unit-level) predictors, and cluster_label is the name of the variable that defines unique clusters. For instance, if w = c("w1","w2), x = "x1", and my clusters are defined by school, then (as far as fitting the model goes) this is equivalent to calling:

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

Problem:

The function runs just fine. My problem is that I want to be able to use anova() to compare two such model objects, but when I try to do this I get the following error, I gather due to my use of reformulate:

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

and, indeed, when I call summary(mod), I get the following line among other returned information:

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

I notice that anova() allows additional arguments to be passed via ..., so is there a way to use this argument so that the function behaves as desired? Or, is there another way to perform a likelihood ratio test on two of my returned model objects without getting this error? Please let me know if I need to provide additional detail.

Reproducible Example:

These are HSB12 data that reproduce the example just fine. This requires the nlme and stringr packages.

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)

Can you please provide a reproducible example to facilitate testing? — Roland, Aug 22 '16 at 14:44
Ah I knew there was a better way to read that .csv file. Thanks for editing! — Jon, Aug 22 '16 at 15:16

score 1 · Accepted Answer · answered Aug 22 '16 at 15:54

There might be a more straight-forward way, but you could use 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.

Using anova() on reformulated() terms in R

1 Answers1