Фон
Я пытаюсь подогнать смешанную модель к функции, основанной на некоторых параметрах. Если я хочу использовать 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
)?
do.call
(хотя, возможно, с теми же побочными эффектами, что и ваше решение). - person aosmith   schedule 17.08.2015