psychometriko psychometriko - 3 months ago 9
R Question

Using anova() on reformulated() terms in R

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)

Answer

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.
Comments