Plot Function For Additive Cox Proportional Hazard Regression


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 Cox-regression 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()
  • non-log 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 y-scale and slight adjustment of the default plotting colors (paler CI shade and stronger term-line)
  • V0.4 – the y-scale 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 (Hat-tip: Arve Ulvik, Eva Pedersen and Roy Nilsen). The option y.log allows both ways (linear and log-scale); 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 x-axis. Any other value for “rug” will omit the rugplot.
  • x.log – logical TRUE/FALSE; log-transformed exposure variable
  • xlab – character; x-axis label
  • ylab – character; y-axis label
  • main – character; main plot title
  • xlim – 2×1 column vector; x-range of plot
  • ylim – 2×1 column vector; y-range of plot
  • col.term – color of HR-curve
  • lwd.term – line width of HR-curve
  • col.se – color of CI (if plotted)
  • cex – numeric; size factor of labels
  • bty – specifies the boxtype around the plot. See ?plot.default
  • Advertisements

14 thoughts on “Plot Function For Additive Cox Proportional Hazard Regression

  1. Thank you for making available this beautiful function plotHR().
    I would like to acknowledge its author(s), who may I mention?
    Vincent

    1. 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 (R-2.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 whois-lookup, my name is Reinhard Seifert, biostatistician at Haukeland University Hospital Bergen, Norway.

      Correspondance: MyForename.MySurname(at)uib.no

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

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

      1. … 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)

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

  3. I really appreciate your contributions! These graphs will likely find a way into my next publication-how shall I reference? Have you considered adding them to bioconductor?

    1. 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 alpha-ish…

      The correct reference for the method I used is:
      Therneau et al. 2000, p. 108
      an excellent reference for COX-regression by the way…

      I just packed the code into a neat function and added some graphical candy, so that my non-geek coworkers can use it…. feel free to use the code šŸ™‚

      Best regards,
      Reinhard

  4. 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 eye-candy 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

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

      1. 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&quot;, 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.

  5. 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 off-topic – I’ve noticed some problems with the prop. hazards assumptions in my data and I want to do a time-spline 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…

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s