Note: This initial blog post is discontinued. If you want to have the latest development of the function have a look at the static page PlotHR
Usually a Coxregression is achieved in R by
library(survival)
model < coxph ( SuvivalObject ~ Covariate1 + Covariate2 + Factor1 + Factor2 , data = Dataset )
The covariates can be enclosed in other funtions:
 factors should be enclosed by factor()
 strata, which allow to adjust for a factor without getting an estimate, should be enclosed by strata()
 nonlog linear continuous terms can be enclosed by
pspline()
In the latter case the model might look like
model < coxph (SurvivalObject ~ pspline(Covariate1) + Covariate2 + factor(Factor1) + strata(Factor2) , data = Dataset )
The functional form of the covariates (including the factors) can now be plotted with
termplot(model)
Though the termplot() function fails with plotting just one covariate and leaves no cusomization.
The function plotHR() plots the functional form of the desired term: plotHR(model)
plots the first term in the model by default but other terms can be accessed by calling their number (e.g. the second one):
plotHR(model , terms = 2)
In order to use the function you have to “source” it into R. It is the same procedure as calling a package, but using “source” instead of “library”.
Paste the function syntax into a textfile and safe it (as plot.HR.R) on your harddisk, remember the path and include
source("C:Path/to/plotHR_0.6.R")
before using the function.
Download plotHR()
Note: I have rewritten the function several times since I wrote the initial post … using version numbers now …

V0.6 – removed the y.log option, since the scale should be logaritmic anyway. Later I will also rewrite the x.log option, since the feature is already incorporated in the plot.default() function. I also removed the dottet line at HR=1 level, since some complained about it overstating the importance of the log(HR) intercept. I included it, since it gives a hint about the significance of the smooth term, in case the confidence intervalls cross over the line… Those who miss it can add manually
lines( h = 0 , type = 2 )
I rewrote the “rugs” option. Try
rugs = "density"
It is still “beta”ish, but some like it.  V0.5 – bug fix for the yscale and slight adjustment of the default plotting colors (paler CI shade and stronger termline)
 V0.4 – the yscale should be logarithmic; a HR of 0.5 (50% reduced Hazard) should show the same distance from HR = 1 as a doubled Hazard (HR = 2); this is now default. The linear scale I used initially is biased in this concern (Hattip: Arve Ulvik, Eva Pedersen and Roy Nilsen). The option y.log allows both ways (linear and logscale); the axis labels denote Hazard Ratio instead of log(HR).
Usage:
plotHR( model , terms = 1 , se = TRUE , rug = "ticks" , x.log = FALSE , xlab = "" , ylab = "Hazard Ratio" , main = NULL , xlim = NULL , ylim = NULL, col.term = "#08519C", lwd.term = 3, col.se = "#DEEBF7", cex = 1 , bty = "n" , axes = TRUE )
 model – a coxph model
 terms – integer; the number of the term to plot
 se – logical TRUE/FALSE; plotting the CI
 rug – “ticks” or “density”; rug plot or density plot at xaxis. Any other value for “rug” will omit the rugplot.
 x.log – logical TRUE/FALSE; logtransformed exposure variable
 xlab – character; xaxis label
 ylab – character; yaxis label
 main – character; main plot title
 xlim – 2×1 column vector; xrange of plot
 ylim – 2×1 column vector; yrange of plot
 col.term – color of HRcurve
 lwd.term – line width of HRcurve
 col.se – color of CI (if plotted)
 cex – numeric; size factor of labels
 bty – specifies the boxtype around the plot. See
?plot.default
Thank you for making available this beautiful function plotHR().
I would like to acknowledge its author(s), who may I mention?
Vincent
Thank you for feedback. I am glad for anybody benefiting from my research work. I really wonder, why there is not anything like it in Therneau’s “survival”package. “termplot()” does not meet the mark… [Update: termplot() has been patched and does work with univariate models (R2.12.1) š ]
Note that I wrote plotHR_0.6.R and it is ready for download.
I would be happy if you reference the blogpost. As you would have seen with a whoislookup, my name is Reinhard Seifert, biostatistician at Haukeland University Hospital Bergen, Norway.
Correspondance: MyForename.MySurname(at)uib.no
I’m having trouble getting this function to work properly. It seems like there is an error extracting the data from the model (to use the words of your comment). As I debug, the value of ‘data’ always returns null. Any thoughts on this?
This looks like an extremely useful function and I’d love to be able to get it to work. Let me know if you think you can help.
Thanks so much!
Hi David,
you have to provide some more detail:
Which version of plotHR do you use?
Did you try
termplot(YourModel , se = T)
to get an impression if everything is alright?
– termplot() is the Goldstandard for plotting the terms. It just does not look nice and fails if there is just one covariate. If termplot() fails, there is something wrong with the model specification.
… you might paste your model specification in a comment which should look like:
library(survival)
model < coxph( YourSurvivalObject ~ pspline(YourSmoothTerm) + SomeOtherCovs , data = YourDataset)
and then of course something like:
source("YourPathTo/plotHR_0.6.R")
plotHR(model)
I discovered my problem was not specifying data= in the coxph statement.
It appears as though this function isn’t compatible with the ns( ) spline, is that correct or am I committing more programming blunders?
Yes, that is true: ns() splines do not work inside coxph(). pspline() only.
I really appreciate your contributions! These graphs will likely find a way into my next publicationhow shall I reference? Have you considered adding them to bioconductor?
Thanks for credits. I have considered suggesting a patch for termplot() or a plot method for coxph() models at some point, but plotHR is quite alphaish…
The correct reference for the method I used is:
Therneau et al. 2000, p. 108
an excellent reference for COXregression by the way…
I just packed the code into a neat function and added some graphical candy, so that my nongeek coworkers can use it…. feel free to use the code š
Best regards,
Reinhard
Thanks for an amazing function. For different reasons I’ve been using the cph in the Design package, I’ve tried to adapt the plotHR to the cph fit but I can’t figure out why I get a:
Error in xy.coords(x, y, xlabel, ylabel, log) :
‘x’ and ‘y’ lengths differ
Right now I’m using coxph for displaying the splines but it’s a little annoying to use different funcitons just for the plot… I you have a good idea to what’s happening it would be awesome if the plotHR was adapted to the cph š
I’ve added some sub() functionality to detect rcs() and other complex formulas. I’ve also added a part for removing interaction terms from the all.labels and some eyecandy for choosing colors for the plot.
I’ve posted my version here: http://pastebin.com/yk0aqUUi
PS I set the pastebin to expire in 30 days
Hi, thanks for contributing. Had a look at your additions and I will definitely consider them.
Concerning your error: I do not even get that far. For me cph() already collides with predictDesign():
Error in predictDesign():
newdata not given and fit did not store x
predictDesign() works with default options, but type = "terms" seems to require newdata to be specified. I just gave it a quick shot today, but even when handing over the dataset as newdata the function complained about different variable lenghts. Looks like missing data related at a glance.
Do you have an idea?
We will figure that one out.
PS.: 'x' and 'y' lengths differ looks very much like missing values on one side of the equation, most of the time one of the covariates missing and plot.default complaining about it…
Another thing: install.packages("Design") threw an error and it turned out that Design is no longer in the CRAN repos. So I dug a bit in and found a forum post (google: “coxph() vs cph()” , sorry no link) of Harrell himself, where he explained that cph() is a wrapper around coxph(), so there seems to be no reason not to use both of them in the same setting.
Regards,
Reinhard
It seems that the “rms” package has succeeded the “Design” package. I got this code to work:
hmohiv<read.table("http://www.ats.ucla.edu/stat/R/examples/asa/hmohiv.csv", sep=",", header = TRUE)
surv < with(hmohiv, Surv(time, censor))
fit < coxph(surv~ pspline(age), data=hmohiv)
plotHR(fit)
ddist < with(hmohiv, datadist(time, age, age, drug))
options(datadist='ddist')
# Passes
fit < cph(surv~ rcs(age, 3), data=hmohiv, x=T, y=T, surv=T)
plotHR(fit, "age")
rcs
# Fails
fit < cph(surv~ rcs(age, 3) + drug, data=hmohiv, x=T, y=T, surv=T)
plotHR(fit, "age")
The problem seems to appear when more variables are added. I'm also a little surprised how different the splines appear between the packages.
Ahh… found the issue: predict for cph ignores the “terms” option. I added a fix for when the term$fit has mutliple column – it may not be beautiful but it works š
I’m a little surprised that the waist of the cph spline is so narrow (se the above example) – I can’t understand quite why this is…
I’ve updated the pastebin to the latest version.
A little offtopic – I’ve noticed some problems with the prop. hazards assumptions in my data and I want to do a timespline the way the coxvc package did but I haven’t had any luck with my post on CrossValidated – perhaps you have an idea? I would like to have a similar plot to the plotHR but for the time axis – the plot that the timecox function provides just seems unappealing to the eye and I believe that the time function in my case is smooth…
Rattling clean web site, appreciate it with this posting.