I've created my own slightly enhanced version of the termplot that I use in this example, you can find it here. I've previously posted on SO but the more I think about it I believe that this probably is more related to the interpretation of the Cox Proportional hazards model than with the actual coding.
The problem
When I look at a Hazard Ratio plot I expect to have a reference point where the confidence interval naturally is 0 and this is the case when I use the cph() from the rms package
but not when I use the coxph() from the survival package
. Is correct behavior by coxph() and if so what is the reference point? Also, the dummy variable in the coxph() has an interval and the value is other than $e^0$?
Example
Here's my test code:
# 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)
The cph plots
This 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")
gives this plot:
The coxph plots
This 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")
gives this plot:
Update
As @Frank Harrell suggested and after adjusting along suggestion in his recent comment I got:
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)))
This gave this very nice plot:
I've looked at the contrast.rms again after the comment and tried this code that gave a plot... although there is probably much more that can be done :-)
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")
Gave this plot:
UPDATE 2
Prof. Thernau was kind enough to comment on the plots lack of a confidence waist:
The smoothing splines in coxph, like the ones in gam, are normalized so that sum(prediction) =0. So I don't have a fixed single point for which the variance is extra small.
Although I'm not yet familiar with GAM this does seem to answer my question: this seems to be an issue of interpretation.
plot
andcontrast
instead ofplot.Predict
andcontrast.rms
. I would useby
orlength
insideseq
instead oftimes
and would givecontrast
two lists so you specify exactly what is being contrasted. You can also use shading withxYplot
for confidence bands. $\endgroup$ – Frank Harrell Feb 6 '12 at 13:06