5

我对我正在尝试编写的函数的行为感到困惑。我的例子来自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包作者 Terry Therneau 的有用意见。这是他的回答:

这是一个源于model.frame做的非标准评估的问题。我发现的唯一方法是将 model.frame=TRUE 添加到原始 coxph 调用中。我认为它是 R 中的一个严重设计缺陷。非标准评估就像是阴暗面——一条诱人而容易的道路,总是以糟糕的结局告终。特里·T。

4

2 回答 2

8

诊断

从错误消息中:

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

请注意,我的 for 示例lm()有点人为,因为我们实际上可以使用newdata参数 inpredict.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论据。文档内容如下:

新数据:

具有与“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这个包是如何工作的。所以我认为我们真的需要一些专家来解释这一点。

于 2016-05-21T15:35:57.003 回答
-2

我想这可能是如果你的

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

作为参数,它没有被正确创建。试试放

is.Surv(formula)

作为函数中的第一行。我怀疑它不起作用,那么我建议使用 apply 系列函数。

于 2016-05-21T17:27:19.923 回答