R - model.frame() und Nicht-Standard-Bewertung

Question

Ich bin verwirrt über ein Verhalten einer Funktion, die ich versuche zu schreiben. Mein Beispiel kommt aus dem survival Paket, aber ich denke, dass die Frage allgemeiner ist. Grundsätzlich ist der Code folgenden

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) 

Wird ein Objekt erhalten, die ich bin interessiert Allerdings, wenn ich es in einer Funktion schreiben,

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) 

die Funktion einen Fehler in der letzten Zeile zurück.:

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) 

Ich bin neugierig ob es etwas offensichtlich ist, dass ich vermisse oder ob solches Verhalten normal ist. Ich finde es seltsam, da ich in der Umgebung von my_function die gleichen Objekte wie in der globalen Umgebung hätte, wenn ich den ersten Teil des Codes ausführe.

Edit: Ich erhielt auch nützliche Informationen von Terry Therneau, dem Autor des survival Pakets. Dies ist seine Antwort:

Dies ist ein Problem, das von der Nicht-Standard-Evaluierung von Model.frame stammt. Der einzige Ausweg, den ich gefunden habe, ist, model.frame = TRUE zum ursprünglichen Coxph-Aufruf hinzuzufügen. Ich halte es für einen ernsthaften Konstruktionsfehler in R. Nicht-Standard-Bewertung ist wie die dunkle Seite - ein verführerischer und einfacher Weg, der immer schlecht endet. Terry T.

Answer 1

Diagnostizieren

Von der Fehlermeldung:

2 survfit(mod_init, newdata = base_case) 
1 my_function(Surv(start, stop, event) ~ rx + number, data = bladder2) 

das Problem eindeutig nicht mit coxph während Modell ist passend, aber mit survfit.

Und von dieser Nachricht:

10 eval(predvars, data, env) 
9 model.frame.default(formula = Surv(start, stop, event) ~ rx + 
    number, data = data) 

Ich kann sagen, dass das Problem ist, dass die Funktion model.frame.default() während der frühen Phase der survfit kann einen Rahmen Modell nicht finden in Formel Surv(start, stop, event) ~ rx + number verwendeten relevante Daten enthält. Daher beklagt es sich.

---

Was ist ein Modellrahmen?

Ein Modellrahmen, aus dem wie lm(), glm() und mgcv:::gam() zu Anpassungsroutine, bestand data Argument gebildet.Es ist ein Datenrahmen mit der gleichen Anzahl von Zeilen wie data, aber:

• fallen alle Variablen nicht durch formula verwiesen
• viele Attribute hinzufügen, die wichtigsten davon envrionement

Die meisten Modellanpassungsroutinen, wie lm(), glm() und mgcv:::gam(), behalten den Modellrahmen standardmäßig in ihrem angepassten Objekt. Dies hat den Vorteil, dass, wenn wir später predict aufrufen und keine newdata bereitgestellt wird, Daten aus diesem Modellrahmen zur Auswertung finden. Ein klarer Nachteil ist jedoch, dass es die Größe Ihres angepassten Objekts wesentlich erhöht.

Allerdings ist survival:::coxph() eine Ausnahme. Es wird standardmäßig nicht solchen Modellrahmen in ihrem angepassten Objekt beibehalten. Nun, das macht das resultierende angepasste Objekt deutlich kleiner, stellt Sie jedoch dem Problem, dem Sie begegnet sind, aus. Wenn wir survival:::coxph() anfordern möchten, um diesen Modellrahmen beizubehalten, dann verwenden Sie model = TRUE dieser Funktion.

---

-Test mit 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) 
    } 

Nun wird dieser Funktionsaufruf fehlschlagen, wie Sie gesehen haben:

my_function(Surv(start, stop, event) ~ rx + number, bladder2, keep.mf = FALSE) 

aber dieser Funktionsaufruf erfolgreich sein wird:

my_function(Surv(start, stop, event) ~ rx + number, bladder2, keep.mf = TRUE) 

---

gleiche Verhalten für lm()

Wir können tatsächlich das gleiche Verhalten in lm() zeigen:

## 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) 
    } 

Jetzt wird dies gelingen, durch Modellrahmen zu halten:

bar(y ~ x - 1, foo, keep.mf = TRUE) 

während dies scheitern wird, b y Verwerfungsmodellrahmen:

bar(y ~ x - 1, foo, keep.mf = FALSE) 

---

Argument Verwendung newdata?

Beachten Sie, dass mein Beispiel für lm() etwas künstlich ist, weil wir eigentlich newdata Argument in predict.lm(), um durch dieses Problem verwenden können:

bar1 <- function(myformula, mydata, keep.mf = TRUE) { 
    fit <- lm(myformula, mydata, model = keep.mf) 
    predict.lm(fit, newdata = lapply(mydata, mean)) 
    } 

nun, ob wir Modellrahmen zu halten, werden beide erfolgreich sein:

bar1(y ~ x - 1, foo, keep.mf = TRUE) 
bar1(y ~ x - 1, foo, keep.mf = FALSE) 

Dann können Sie sich fragen: Können wir das gleiche für survfit() tun?

survfit() ist eine generische Funktion, in Ihrem Code rufen Sie wirklich survfit.coxph(). Es gibt tatsächlich ein newdata Argument für diese Funktion. Die Dokumentation lautet:

newdata:

ein Datenrahmen mit denselben Variablennamen wie diejenigen, die in der ‚coxph‘ Formel erscheinen. ... ... Default ist der Mittelwert der Kovariaten, die in der 'Coxph'-Anpassung verwendet werden.

Also, lassen Sie uns versuchen:

my_function1 <- function(myformula, mydata) { 
    mtrace.off() 
    fit <- coxph(myformula, mydata, method = "breslow") 
    survival:::survfit.coxph(fit, newdata = lapply(mydata, mean)) 
    } 

und wir hoffen, dass diese Arbeit:

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

Aber:

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) 

Beachten Sie, dass, obwohl wir in newdata passieren, ist es nicht im Bau von Modellrahmen verwendet:

3: stats::model.frame(object) 

Nur object, eine Kopie des angepassten Modells wird auf model.frame.default() weitergegeben.

Dies ist sehr unterschiedlich zu dem, was in predict.lm(), predict.glm() und mgcv:::predict.gam() passiert. In diesen Routinen wird newdata an model.frame.default() übergeben. Zum Beispiel in lm() gibt es:

m <- model.frame(Terms, newdata, na.action = na.action, xlev = object$xlevels) 

Ich benutze kein survival Paket, so dass nicht sicher, wie newdata Arbeiten in diesem Paket. Ich denke, wir brauchen wirklich einen Experten, der das erklärt.

Answer 2

ich denke, es könnte sein, dass, wenn Ihr

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

als Parameter ist, ist es nicht richtig erstellt wird, erhalten. Versuchen Sie setzen

is.Surv(formula) 

als Ihre erste Zeile in der Funktion. Ich vermute, dass es nicht funktioniert, dann würde ich empfehlen, eine Familie von Funktionen zu verwenden.