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)

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

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)

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