# Different prediction plot from survival coxph and rms cph

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

-
 One for @Frank Harrell ? – onestop Jan 30 '12 at 19:35 Several comments. First read biostat.mc.vanderbilt.edu/Rrms for differences between rms and Design packages. Second, use plot() instead of plot.Predict to save work. Third you can easily generate plots for both sexes, e.g. using Predict(fit, age, sex, fun=exp) # exp=anti-log; then plot(result) or plot(result, ~ age | sex). You don't use "x=NA" in Predict. rms uses lattice graphics so usual par graphics parameters and mfrow don't apply. See examples in my rms course handout at biostat.mc.vanderbilt.edu/rms. For contrast.rms study the documentation more. – Frank Harrell Jan 31 '12 at 22:24 Thank you very much for your input. I've updated the code with better examples and added prof. Thernau's response. PS I'm really excited that your planning a new version of the book, extending the cutpoint bias section will be very useful as a reference – Max Gordon Feb 6 '12 at 5:04 You can use plot and contrast instead of plot.Predict and contrast.rms. I would use by or length inside seq instead of times and would give contrast two lists so you specify exactly what is being contrasted. You can also use shading with xYplot for confidence bands. – Frank Harrell Feb 6 '12 at 13:06 Thanks. I like using the plot.Predict because then I get to the right help in RStudio - something that in my case is much more vital than the time it takes to write the full function name (by using autocomplete (tab) I actually don't loose that much time). – Max Gordon Feb 6 '12 at 16:44