Different prediction plot from survival coxph and rms cph

Question

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.

Answer 1

I think there should definitely be a point where the confidence interval is zero width. You might also try a third way which is to use solely rms functions. There is an example under the help file for contrast.rms to get a hazard ratio plot. It starts with the comment # show separate estimates by treatment and sex. You'll need to anti-log to get the ratio.