R - model.frame() и нестандартная оценка

Question

Я озадачен поведением функции, которую я пытаюсь написать. Мой пример исходит из пакета survival, но я думаю, что вопрос более общий. В основном, следующий код

library(survival) 
data(bladder) ## this will load "bladder", "bladder1" and "bladder2" 

mod_init <- coxph(Surv(start, stop, event) ~ rx + number, data = bladder2, method = "breslow") 
survfit(mod_init) 

даст объект, который я заинтересован в Тем не менее, когда я пишу это в функции,

my_function <- function(formula, data) { 
    mod_init <- coxph(formula = formula, data = data, method = "breslow") 
    survfit(mod_init) 
    } 

my_function(Surv(start, stop, event) ~ rx + number, data = bladder2) 

функция возвращает ошибку в последней строке.:

Error in eval(predvars, data, env) : 
    invalid 'envir' argument of type 'closure' 
10 eval(predvars, data, env) 
9 model.frame.default(formula = Surv(start, stop, event) ~ rx + 
    number, data = data) 
8 stats::model.frame(formula = Surv(start, stop, event) ~ rx + 
    number, data = data) 
7 eval(expr, envir, enclos) 
6 eval(temp, environment(formula$terms), parent.frame()) 
5 model.frame.coxph(object) 
4 stats::model.frame(object) 
3 survfit.coxph(mod_init) 
2 survfit(mod_init) 
1 my_function(Surv(start, stop, event) ~ rx + number, data = bladder2) 

Мне любопытно, есть ли что-то очевидное, что мне не хватает или такое поведение нормальное. Я нахожу это странным, поскольку в среде my_function у меня были бы те же объекты, что и в глобальной среде при запуске первой части кода.

Редактировать: Я также получил полезный вклад от Терри Терно, автора пакета survival. Это его ответ:

Это проблема, связанная с нестандартной оценкой, выполненной с помощью model.frame. Единственный выход из этого, что я нашел, это добавить model.frame = TRUE в исходный вызов coxph. Я считаю это серьезным недостатком дизайна в R. Нестандартная оценка похожа на темную сторону - заманчивый и легкий путь, который всегда заканчивается плохо. Терри Т.

Answer 1

Диагностируйте

С сообщением об ошибке:

2 survfit(mod_init, newdata = base_case) 
1 my_function(Surv(start, stop, event) ~ rx + number, data = bladder2) 

проблема явно не с coxph во время подгонки модели, но с survfit.

И из этого сообщения:

10 eval(predvars, data, env) 
9 model.frame.default(formula = Surv(start, stop, event) ~ rx + 
    number, data = data) 

Я могу сказать, что проблема заключается в том, что во время ранней стадии survfit, функция model.frame.default() не может найти модели кадра, содержащий соответствующие данные, используемые в формуле Surv(start, stop, event) ~ rx + number. Поэтому он жалуется.

---

Что такое рама модели?

модели кадр, формируются из data аргумента, переданного фитинг рутины, как lm(), glm() и mgcv:::gam().Это кадр данных с тем же числом строк, как data, но:

• сбросив все переменные не ссылаются formula
• добавив множество атрибутов, наиболее важным из которых является envrionement

Наиболее подпрограммы модели, такие как lm(), glm() и mgcv:::gam(), по умолчанию сохраняют рамку модели в соответствующем объекте. Это имеет преимущество, что, если мы позже позвоним predict, и не будет newdata, он найдет данные из этого модельного кадра для оценки. Однако явным недостатком является то, что он существенно увеличит размер вашего приспособленного объекта.

Однако, исключение составляет survival:::coxph(). Он по умолчанию не сохранит такую ​​модельную рамку в установленном объекте. Ясно, что это делает результирующий приспособленный объект намного меньшим по размеру, но подвергает вас проблеме, с которой вы столкнулись. Если мы хотим задать survival:::coxph(), чтобы сохранить эту модель, используйте model = TRUE этой функции.

---

Тест с survial:::coxph()

library(survival); data(bladder) 

my_function <- function(myformula, mydata, keep.mf = TRUE) { 
    fit <- coxph(myformula, mydata, method = "breslow", model = keep.mf) 
    survfit(fit) 
    } 

Теперь этот вызов функции не получится, как вы уже видели:

my_function(Surv(start, stop, event) ~ rx + number, bladder2, keep.mf = FALSE) 

, но этот вызов функции будет иметь успех:

my_function(Surv(start, stop, event) ~ rx + number, bladder2, keep.mf = TRUE) 

---

То же поведение для lm()

Мы можем на самом деле демонстрируют такое же поведение в lm():

## generate some toy data 
foo <- data.frame(x = seq(0, 1, length = 20), y = seq(0, 1, length = 20) + rnorm(20, 0, 0.15)) 

## a wrapper function 
bar <- function(myformula, mydata, keep.mf = TRUE) { 
    fit <- lm(myformula, mydata, model = keep.mf) 
    predict.lm(fit) 
    } 

Теперь это будет успех, сохраняя модель рамы:

bar(y ~ x - 1, foo, keep.mf = TRUE) 

пока это не удастся, b у отбрасывания модель рамы:

bar(y ~ x - 1, foo, keep.mf = FALSE) 

---

Использование аргумента newdata?

Обратите внимание, что мой пример для lm() немного искусственно, потому что мы можем фактически использовать newdata аргумент в predict.lm(), чтобы пройти через эту проблему:

bar1 <- function(myformula, mydata, keep.mf = TRUE) { 
    fit <- lm(myformula, mydata, model = keep.mf) 
    predict.lm(fit, newdata = lapply(mydata, mean)) 
    } 

Теперь держите ли мы модель кадр, как получится:

bar1(y ~ x - 1, foo, keep.mf = TRUE) 
bar1(y ~ x - 1, foo, keep.mf = FALSE) 

Тогда вы можете задаться вопросом: можем ли мы сделать то же самое для survfit()?

survfit() - общая функция, в вашем коде вы действительно звоните survfit.coxph(). Для этой функции действительно есть аргумент newdata. Документация гласит:

NewData:

фрейм данных с одинаковыми именами переменных, как те, которые появляются в «coxph» формулы . ... ... Значение по умолчанию - среднее значение ковариатов, используемых в соответствии с 'coxph'.

Итак, давайте попробуем:

my_function1 <- function(myformula, mydata) { 
    mtrace.off() 
    fit <- coxph(myformula, mydata, method = "breslow") 
    survival:::survfit.coxph(fit, newdata = lapply(mydata, mean)) 
    } 

, и мы надеемся, что эта работа:

my_function1(Surv(start, stop, event) ~ rx + number, bladder2) 

Но:

Error in is.data.frame(data) (from #5) : object 'mydata' not found 

1: my_function1(Surv(start, stop, event) ~ rx + number, bladder2) 
2: #5: survival:::survfit.coxph(fit, lapply(mydata, mean)) 
3: stats::model.frame(object) 
4: model.frame.coxph(object) 
5: eval(temp, environment(formula$terms), parent.frame()) 
6: eval(expr, envir, enclos) 
7: stats::model.frame(formula = Surv(start, stop, event) ~ rx + number, data = 
8: model.frame.default(formula = Surv(start, stop, event) ~ rx + number, data 
9: is.data.frame(data) 

Обратите внимание, что хотя мы проходим в newdata, то есть не используется при построении рамы модели:

3: stats::model.frame(object) 

Только object, копия подобранной модели, передается model.frame.default().

Это очень отличается от того, что происходит в predict.lm(), predict.glm() и mgcv:::predict.gam(). В этих подпрограммах newdata передается model.frame.default(). Например, в lm(), есть:

m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels) 

Я не использую survival пакет, так что не знаю, как newdata работы в этом пакете. Поэтому я думаю, что нам действительно нужен эксперт, объясняющий это.

Answer 2

Я думаю, что это может быть, что если ваш

Surv(start, stop, event) ~ rx + number 

в качестве параметра, он не получает правильно создан. Попробуйте поместить

is.Surv(formula) 

как ваша первая линия в функции. Я подозреваю, что это не сработает, тогда я бы предложил использовать семейство функций.