12
$\begingroup$

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:

cph() termplot2

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:

coxph() termplot2

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:

Lattice 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:

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

Improve this question
$\endgroup$
4
  • 5
    $\begingroup$ 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. $\endgroup$ Commented Jan 31, 2012 at 22:24
  • 1
    $\begingroup$ 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 $\endgroup$ Commented Feb 6, 2012 at 5:04
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 6, 2012 at 13:06
  • 1
    $\begingroup$ 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). $\endgroup$ Commented Feb 6, 2012 at 16:44

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thank you for your answer. Do you think I should mention this problem to prof. Terry Therneau if it is to be considered as a bug/misinterpretation? I've also looked into the graphic solutions in the rms package, I can't quite understand the usage of contrast.rms for plots. The plot.Predict seems to do a termplot similar output but I can't get it to do exactly what I want... see my update to the question. $\endgroup$ Commented Jan 31, 2012 at 21:35
  • 5
    $\begingroup$ It would be good to write him to enquire, and tell him thanks for the ride to the airport he gave me a few minutes ago. I'll comment above on the other questions. $\endgroup$ Commented Jan 31, 2012 at 22:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.