Unterschiedliches Vorhersageplot von Survival Coxph und RMS Cph

9

Ich habe meine eigene, leicht verbesserte Version des Termplots erstellt, die ich in diesem Beispiel verwende. Sie finden sie hier . Ich habe zuvor auf SO gepostet, aber je mehr ich darüber nachdenke, desto mehr hängt dies wahrscheinlich mit der Interpretation des Cox Proportional Hazards-Modells zusammen als mit der tatsächlichen Codierung.

Das Problem

Wenn ich mir ein Hazard Ratio-Diagramm anschaue, erwarte ich einen Referenzpunkt, an dem das Konfidenzintervall natürlich 0 ist, und dies ist der Fall, wenn ich cph () aus dem verwende, rms packageaber nicht, wenn ich coxph () aus dem verwende survival package. Ist korrektes Verhalten von coxph () und wenn ja, was ist der Bezugspunkt? Außerdem hat die Dummy-Variable in coxph () ein Intervall und der Wert ist anders als ?e0

Beispiel

Hier ist mein Testcode:

# Load libs
library(survival)
library(rms)

# Regular survival
survobj <- with(lung, Surv(time,status))

# Prepare the variables
lung$sex <- factor(lung$sex, levels=1:2, labels=c("Male", "Female"))
labels(lung$sex) <- "Sex"
labels(lung$age) <- "Age"

# The rms survival
ddist <- datadist(lung)
options(datadist="ddist")
rms_surv_fit <- cph(survobj~rcs(age, 4)+sex, data=lung, x=T, y=T)

Die cph-Diagramme

Dieser Code:

termplot2(rms_surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05,
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("cph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

gibt diese Handlung:

cph () termplot2

Die Coxph-Diagramme

Dieser Code:

termplot2(surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05, 
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("coxph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

gibt diese Handlung:

coxph () termplot2

Aktualisieren

Wie @Frank Harrell vorschlug und nachdem ich den Vorschlag in seinem letzten Kommentar angepasst hatte, bekam ich:

p <- Predict(rms_surv_fit, age=seq(50, 70, times=20), 
             sex=c("Male", "Female"), fun=exp)
plot.Predict(p, ~ age | sex,
             col="black",
             col.fill=gray(seq(.8, .75, length=5)))

Dies gab diese sehr schöne Handlung:

Gitterplot

Ich habe mir nach dem Kommentar noch einmal die kontrast.rms angesehen und diesen Code ausprobiert, der eine Handlung ergab ... obwohl wahrscheinlich noch viel mehr getan werden kann :-)

w <- contrast.rms(rms_surv_fit, 
                  list(sex=c("Male", "Female"), 
                       age=seq(50, 70, times=20)))

xYplot(Cbind(Contrast, Lower, Upper) ~ age | sex, 
       data=w, method="bands")

Gab diese Handlung:

Das Kontrastdiagramm

UPDATE 2

Prof. Thernau war so freundlich, die Handlungen zu kommentieren, denen es an Selbstvertrauen mangelt:

Die Glättungssplines in Coxph werden wie die in Gam so normalisiert, dass die Summe (Vorhersage) = 0 ist. Ich habe also keinen festen Einzelpunkt, für den die Varianz besonders klein ist.

Obwohl ich mit GAM noch nicht vertraut bin, scheint dies meine Frage zu beantworten: Dies scheint eine Frage der Interpretation zu sein.

Max Gordon
quelle
3
Mehrere Kommentare. Lesen Sie zuerst biostat.mc.vanderbilt.edu/Rrms, um Unterschiede zwischen RMS- und Design-Paketen zu ermitteln. Zweitens verwenden Sie plot () anstelle von plot.Predict, um Arbeit zu speichern. Drittens können Sie auf einfache Weise Diagramme für beide Geschlechter erstellen, z. B. mithilfe von Predict (Passform, Alter, Geschlecht, Spaß = exp) # exp = anti-log; dann Handlung (Ergebnis) oder Handlung (Ergebnis, ~ Alter | Geschlecht). Sie verwenden "x = NA" in Predict nicht. rms verwendet Gittergrafiken, sodass die üblichen Par-Grafikparameter und mfrow nicht zutreffen. Beispiele finden Sie in meinem Handout zum Effektivkurs unter biostat.mc.vanderbilt.edu/rms . Für Kontrast.rms studieren Sie die Dokumentation mehr.
Frank Harrell
1
Vielen Dank für Ihre Eingabe. Ich habe den Code mit besseren Beispielen aktualisiert und prof hinzugefügt. Thernaus Antwort. PS Ich freue mich sehr, dass Ihre Planung einer neuen Version des Buches, die den Abschnitt mit den Cutpoint-Bias erweitert, als Referenz sehr nützlich sein wird
Max Gordon
1
Sie können plotund contrastanstelle von plot.Predictund verwenden contrast.rms. Ich würde byoder lengthinnen seqanstelle von timesund würde contrastzwei Listen geben, damit Sie genau angeben, was kontrastiert wird. Sie können auch Schattierungen mit xYplotfür Vertrauensbereiche verwenden.
Frank Harrell
1
Vielen Dank. Ich benutze gerne den Plot.Predict, weil ich dann in RStudio die richtige Hilfe bekomme - etwas, das in meinem Fall viel wichtiger ist als die Zeit, die zum Schreiben des vollständigen Funktionsnamens benötigt wird (mithilfe der automatischen Vervollständigung (Registerkarte), die ich eigentlich nicht habe) so viel Zeit verlieren).
Max Gordon

Antworten:

5

Ich denke, es sollte definitiv einen Punkt geben, an dem das Konfidenzintervall die Breite Null hat. Sie können auch einen dritten Weg ausprobieren, bei dem ausschließlich Effektivfunktionen verwendet werden. Unter der Hilfedatei befindet sich ein Beispiel für kontrast.rms, um ein Hazard-Ratio-Diagramm zu erhalten. Es beginnt mit dem Kommentar # zeigt separate Schätzungen nach Behandlung und Geschlecht. Sie müssen Anti-Log, um das Verhältnis zu erhalten.

Frank Harrell
quelle
1
Vielen Dank für Ihre Antwort. Glaubst du, ich sollte dieses Problem prof. Terry Therneau, wenn es als Fehler / Fehlinterpretation zu betrachten ist? Ich habe mir auch die Grafiklösungen im RMS-Paket angesehen. Ich kann die Verwendung von Kontrast.RMS für Diagramme nicht ganz verstehen. The plot.Predict scheint eine Termplot-ähnliche Ausgabe zu machen, aber ich kann es nicht dazu bringen, genau das zu tun, was ich will ... siehe mein Update der Frage.
Max Gordon
2
Es wäre gut, ihm zu schreiben, um sich zu erkundigen und ihm für die Fahrt zum Flughafen zu danken, die er mir vor ein paar Minuten gegeben hat. Ich werde oben auf die anderen Fragen eingehen.
Frank Harrell