診斷
從錯誤消息：
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'fit中使用的協變量的均值。
所以，讓我們試試：
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作品在這個包怎麼樣。所以我認爲我們確實需要一些專家解釋這一點。
謝謝你非常清楚的解釋。 但有一件事是困擾着我：爲什麼在交互式使用中工作（在函數之外）？ – Theodor