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 package
mais 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 ?
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:
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:
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:
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:
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.
plot
etcontrast
au lieu deplot.Predict
etcontrast.rms
. J'utiliseraisby
ou à l'length
intérieurseq
au lieu detimes
et donneraiscontrast
deux listes afin que vous spécifiez exactement ce qui est contrasté. Vous pouvez également utiliser l'ombrage avecxYplot
des bandes de confiance.Réponses:
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.
la source