Создать объект 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, чтобы их можно было получить также в разных средах.

Я могу жить со своим обходным путем, хотя он имеет некоторые уродливые побочные эффекты, поскольку print(mm1) показывает все данные вместо простых имен. Итак, мой вопрос в том, есть ли какие-то другие стратегии, чтобы получить то, что я намереваюсь? Или я должен написать сопровождающему contrast и спросить его, может ли он изменить код для объектов lme так, чтобы он больше не полагался на слот call, а пытался решить проблему иначе (как это делается для lm)?


person thothal    schedule 03.08.2015    source источник
comment
Я думаю, что это может быть дубликат this question, с ответом, показывающим, как это сделать через do.call (хотя, возможно, с теми же побочными эффектами, что и ваше решение).   -  person aosmith    schedule 17.08.2015
comment
Да, я видел это. Проблема не столько в том, что я не могу ее решить, вопрос был в том, как избежать побочных эффектов.   -  person thothal    schedule 18.08.2015


Ответы (1)


Я думаю, что то, с чем вы боретесь, - это просто ошибочная реализация 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
person Forrest R. Stevens    schedule 18.08.2015
comment
Спасибо, это действительно помогает, и я также узнал кое-что новое (assignInNamespace). Более чем стоит награда :) Спасибо. - person thothal; 19.08.2015
comment
Хех, рад помочь. Это маловероятно в вашем случае, но возня с функциями внутри пространств имен пакетов может привести к некоторым непредвиденным последствиям (например, если по какой-то причине contrastCalc() повторно вызывает contrast.lme()). Но иногда это очень полезно, особенно если вы тестируете незначительные изменения и не хотите иметь дело с исходным кодом всего пакета. - person Forrest R. Stevens; 20.08.2015