诊断
从错误消息中:
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
这个包是如何工作的。所以我认为我们真的需要一些专家来解释这一点。