Create lme object inside function
Background
I am trying to set 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, since it contrast
uses the slot call
from the object lme
to define the arguments data
, fixed
or random
is passed to lme
in functions (see code). This, by the way, does not apply to the object lm
.
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 traced the error down to the part where contrast
evaluating the slot fixed
in the slot call
mm2
is equal mF
, which of course is not known at the top level since it was only defined inside my function makeMixedModel2
. The workaround in this makeMixedModel1
eliminates the fact that by explicitly rewriting the appropriate slots in call
.
Apparently, for objects lm
this is solved in a more reasonable way, since there is no need to perform manual overwrite, since contrast
it seems to tempt all the parts in the correct context, although, of course, mF
and mdat
also not known:
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 my guess is that it lm
stores values formula
and data
somewhere, so in can be retrieved in different environments as well.
I can live with my workaround, although it has some ugly side effects as it print(mm1)
shows all data instead of simple names. So my question is, are there any other strategies for getting what I intend? Or do I need to write to a co-worker contrast
and ask him if he can change the code for the objects lme
so that hes no longer rely on the slot call
but tries to fix the problem otherwise (how is it done for lm
)?
source to share
I think what you are fighting is just a mistake contrast()
for lme
objects. I would contact the author to fix it (this may be the result of a recent change with help nlme
). But at the same time, you can avoid side effects by implementing your workaround in a function contrast.lme()
rather than a model building 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))
Productivity:
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)
Productivity:
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
source to share