r - 在函数内创建lme对象

背景

我正在尝试根据一些参数将混合模型拟合到函数中。如果要使用contrast中的library(contrast)，则必须使用一种解决方法，因为contrast使用call对象中的lme插槽来确定data，fixed或random参数传递给函数中的lme（参见代码）。顺便说一下，lm对象不是这种情况。

数据

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

码

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

题

我已将错误跟踪到contrast评估fixed的call插槽内的mm2插槽的部分，该部分等于mF当然是在顶层，因为它已定义仅在我的函数makeMixedModel2中。 makeMixedModel1中的解决方法通过显式覆盖call中的相应插槽来进行补救。

显然，对于lm对象，这可以通过更智能的方式解决，因为无需进行手动覆盖，因为contrast似乎在正确的上下文中撤回了所有部分，尽管mF和mdat当然不是已知：

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))

因此，我假设lm将formula和data的值存储在某个位置，以便可以在不同的环境中检索in。

我可以忍受我的解决方法，尽管它有一些丑陋的副作用，因为print(mm1)显示所有数据而不是简单的名称。所以我的问题是，是否有其他一些策略可以达到我的预期？还是我必须写信给contrast的维护者并问他，他是否可以更改lme对象的代码，使他不再依赖call插槽，而是尝试解决该问题（如对lm所做的那样？

最佳答案

我认为您要解决的只是contrast()对象的lme错误实现。我会联系作者进行修复（这可能是最近使用nlme进行更改的结果）。但是与此同时，您可以通过在contrast.lme()函数而不是模型构造函数中实现变通办法来避免副作用：

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))

产量：

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

和：

print(mm2)

产量：

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

关于r - 在函数内创建lme对象，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/31788607/