诊断
来自错误信息：
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() 示例有点人为，因为我们实际上可以在predict.lm() 中使用newdata 参数来解决这个问题：
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 在这个包中是如何工作的。所以我认为我们真的需要一些专家来解释这一点。