Create lme object within a function

Question

Background

I am trying to fit a mixed model in a function based on some parameters. If I want to use contrast from library(contrast) I have to use a workaround, as contrast uses the call slot from the lme object to determine the data, fixed or random argument passed to lme in the function (cf. code). This is by the way not the case with an lm object.

Data

set.seed(1)
dat <- data.frame(x = runif(100), y1 = rnorm(100), y2 = rnorm(100),
                  grp = factor(sample(2, 100, replace = TRUE)))

Code

library(contrast)
library(nlme)
makeMixedModel1 <- function(resp, mdat = dat) {
   mF <- formula(paste(resp, "x", sep = "~"))
   mdat <- mdat[resp > 0, ]
   mod <- lme(mF, random = ~ 1 | grp, data = mdat)
   mC <- mod$call
   mC$fixed <- mF
   mC$data <- mdat
   mod$call <- mC
   mod
}

makeMixedModel2 <- function(resp, mdat = dat) {
   mF <- formula(paste(resp, "x", sep = "~"))
   mdat <- mdat[resp > 0, ]
   lme(mF, random = ~ 1 | grp, data = mdat)
}

mm1 <- makeMixedModel1("y1")
mm2 <- makeMixedModel2("y1")
contrast(mm1, list(x = 1)) ## works as expected
# lme model parameter contrast
# 
#   Contrast      S.E.      Lower     Upper    t df Pr(>|t|)
#  0.1613734 0.2169281 -0.2692255 0.5919722 0.74 96   0.4588

contrast(mm2, list(x = 1)) ## gives an error
# Error in eval(expr, envir, enclos) : object 'mF' not found

Question

I have tracked down the error to the part where contrast evaluates the fixedslot within the call slot of mm2 whis equals to mF which is of course not known at the top level, as it was defined only inside my function makeMixedModel2. The workaround in makeMixedModel1 remedies that by explicitely overwriting the respective slots in call.

Apparently, for lm objects this is solved in a smarter way, as there's no need to do the manual overwriting as contrast seems to evaulate all the parts within the right context though of course mF and mdat are not known either:

makeLinearModel <- function(resp, mdat = dat) {
   mF <- formula(paste(resp, "x", sep = "~"))
   mdat <- mdat[resp > 0, ]
   lm(mF, data = mdat)
}
contrast(makeLinearModel("y1"), list(x = 1))

So, I assume that lm stores the values of the formula and the data somewhere, such that in can be retrieved also at different environments.

I can live with my workaround, though it has some ugly side-effects as a print(mm1) shows all the data instead of a simple names. So my question is, whether there are some other strategies to get what I intend? Or do I have to write to the maintainer of contrast and ask him, whether he can change the code for lme objects such that hes does not rely anymore on the call slot, but tries to solve the issue otherwise (as it is done for the lm)?

Answer 1

I think what you're fighting is just a buggy implementation of contrast() for lme objects. I would contact the author to fix it (it may be a result of something recent changing with nlme). But in the mean time you can avoid your side effects by implementing your workaround in the contrast.lme() function rather than in your model construction function:

contrast.lme <- function(fit, ...) {
   mC <- fit$call
   mC$fixed <- formula(fit) 
   mC$data <- fit$data
   fit$call <- mC

   library(nlme)
   contrast:::contrastCalc(fit, ...)
}
assignInNamespace("contrast.lme", contrast.lme, "contrast")

mm2 <- makeMixedModel2("y1")

contrast(mm2, list(x = 1))

Yields:

lme model parameter contrast

  Contrast      S.E.      Lower     Upper    t df Pr(>|t|)
 0.1613734 0.2169281 -0.2692255 0.5919722 0.74 96   0.4588

And:

print(mm2)

Yields:

Linear mixed-effects model fit by REML
  Data: mdat 
  Log-restricted-likelihood: -136.2472
  Fixed: mF 
(Intercept)           x 
 -0.1936347   0.3550081 

Random effects:
 Formula: ~1 | grp
        (Intercept)  Residual
StdDev:    0.131666 0.9365614

Number of Observations: 100
Number of Groups: 2