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.
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$?
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
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
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:
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:
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.