Graphique de prédiction différent de coxph de survie et rms cph

9

J'ai créé ma propre version légèrement améliorée du termplot que j'utilise dans cet exemple, vous pouvez le trouver ici . J'ai déjà posté sur SO mais plus j'y pense, plus je pense que cela est probablement plus lié à l'interprétation du modèle des risques proportionnels de Cox qu'au codage réel.

Le problème

Quand je regarde un tracé de Hazard Ratio, je m'attends à avoir un point de référence où l'intervalle de confiance est naturellement 0 et c'est le cas lorsque j'utilise le cph () du rms packagemais pas quand j'utilise le coxph () du survival package. Le comportement est-il correct par coxph () et si oui quel est le point de référence? En outre, la variable fictive dans le coxph () a un intervalle et la valeur est différente de ?e0

Exemple

Voici mon code de test:

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

Les parcelles cph

Ce 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")

donne ce complot:

cph () termplot2

Les intrigues coxph

Ce 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")

donne ce complot:

coxph () termplot2

Mise à jour

Comme l'a suggéré @Frank Harrell et après avoir ajusté la suggestion dans son récent commentaire, j'ai eu:

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

Cela a donné cette très belle intrigue:

Tracé en treillis

J'ai de nouveau regardé contrast.rms après le commentaire et essayé ce code qui a donné un tracé ... bien qu'il y ait probablement beaucoup plus à faire :-)

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")

A donné ce complot:

L'intrigue de contraste

MISE À JOUR 2

Le professeur Thernau a eu la gentillesse de commenter le manque de confiance des complots:

Les splines de lissage en coxph, comme celles de gam, sont normalisées de sorte que la somme (prédiction) = 0. Je n'ai donc pas de point unique fixe pour lequel la variance est très petite.

Bien que je ne sois pas encore familier avec GAM, cela semble répondre à ma question: cela semble être un problème d'interprétation.

Max Gordon
la source
3
Plusieurs commentaires. Lisez d'abord biostat.mc.vanderbilt.edu/Rrms pour les différences entre les packages rms et Design. Deuxièmement, utilisez plot () au lieu de plot.Predict pour économiser du travail. Troisièmement, vous pouvez facilement générer des graphiques pour les deux sexes, par exemple en utilisant Predict (fit, age, sex, fun = exp) # exp = anti-log; puis tracer (résultat) ou tracer (résultat, ~ âge | sexe). Vous n'utilisez pas "x = NA" dans Predict. rms utilise des graphiques en treillis, les paramètres graphiques par et mfrow ne s'appliquent pas. Voir des exemples dans mon document de cours rms à biostat.mc.vanderbilt.edu/rms . Pour contrast.rms, étudiez davantage la documentation.
Frank Harrell
1
Merci beaucoup pour votre contribution. J'ai mis à jour le code avec de meilleurs exemples et ajouté prof. Réponse de Thernau. PS Je suis vraiment excité que votre planification d'une nouvelle version du livre, étendant la section du biais de point de coupure soit très utile comme référence
Max Gordon
1
Vous pouvez utiliser plotet contrastau lieu de plot.Predictet contrast.rms. J'utiliserais byou à l' lengthintérieur seqau lieu de timeset donnerais contrastdeux listes afin que vous spécifiez exactement ce qui est contrasté. Vous pouvez également utiliser l'ombrage avec xYplotdes bandes de confiance.
Frank Harrell
1
Merci. J'aime utiliser l'intrigue.Prédict parce qu'alors j'obtiens l'aide appropriée dans RStudio - quelque chose qui dans mon cas est beaucoup plus vital que le temps qu'il faut pour écrire le nom complet de la fonction (en utilisant la saisie semi-automatique (tab), je ne le fais pas vraiment perdre autant de temps).
Max Gordon

Réponses:

5

Je pense qu'il devrait certainement y avoir un point où l'intervalle de confiance est de largeur nulle. Vous pouvez également essayer une troisième méthode qui consiste à utiliser uniquement des fonctions RMS. Il y a un exemple sous le fichier d'aide pour contrast.rms pour obtenir un tracé du rapport de risque. Cela commence par le commentaire # montrer des estimations distinctes par traitement et par sexe. Vous aurez besoin d'anti-log pour obtenir le ratio.

Frank Harrell
la source
1
Merci pour votre réponse. Pensez-vous que je devrais mentionner ce problème au prof. Terry Therneau si cela doit être considéré comme un bug / une mauvaise interprétation? J'ai également examiné les solutions graphiques du package rms, je ne comprends pas très bien l'utilisation de contrast.rms pour les tracés. L'intrigue.Predict semble faire une sortie similaire à l'intrigue, mais je n'arrive pas à faire exactement ce que je veux ... voir ma mise à jour de la question.
Max Gordon
2
Ce serait bien de lui écrire pour s'enquérir et lui dire merci pour le trajet jusqu'à l'aéroport qu'il m'a donné il y a quelques minutes. Je commenterai ci-dessus les autres questions.
Frank Harrell