cox

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
    • Advertisement

Additive COX-regression

Update:I have written a much more detailed static page about the additive COX model: http://rforge.org/plothr/
The page has a download link to the function plotHR() which does all the fuzz. It is extensively commented. It should be easy to understand the syntax and modify it for individual purposes.

Therneau et al. refer to the proportional hazards model or COX-regression model as “the workhorse of regression analysis for censored data”. They show how to implement the additive form of this model in SAS and S-pluss; already mentioned by Hastie and Tibshirany in 1986 when introducing Generalized Additive Models (GAM).

I found modelling the functional form of the covariates in a regression model for rightcensored survival times with smoothing splines extremely useful. And the implementation is absolutely straightforward in R.

The only thing needed is the installation of the R-libraries “survival” and “pspline”:

install.packages("pspline")
and
install.packages("survival")

In the following code I will refer to a dataset “MyData” with a binary status variable “death” and a time-to-event variable “days2death”.
The status variable “death” should be (not necessarily) 1 if the event of interesst occured to the subject and “days2death” gives then the time to this event.

Viualizing the functional form of a covariate takes the following steps:

  1. create the survival object of interesst
  2. fit a proportional hazards model with smoothing splines,
  3. predict the functional form of the covariate of interesst and
  4. plot it!

Note that there is the termplot() function in R which gives you the GAM plots after the modelfit, so step 3 would not be necessary – BUT: it has a bug and fails plotting a single covariate; and it does not allow all to much customizing.

This is the R code to achieve the analysis:

1 Create survival object:

surv.death <- Surv(MyData$days2death, MyData$death)

2 Fit proportional hazards model with smoothing splines for continuous covariates:

library(survival)
library(pspline)
pham.fit <- coxph( surv.death ~ pspline(EF, df=4) + pspline(Age, df=4) + strata (Sex, df=4) , data = MyData)

The model above includes the continuous covariates “EF” (ejection fraction) and “Age” and stratifies for “Sex”.

3 Produce the fitted smoothing spline for the first covariate in the above model formula with standard errors

predicted <- predict(pham.fit , type = "terms" , se.fit = TRUE , terms = 1)
“terms=1” refers to “pspline(EF,df=4)”

4 Plot it

First plotting axes and labels
plot(0 , xlab="Ejection Fraction" , ylab = "Hazard Ratio" , main = "All-cause Death" , type = "n" , xlim=c(0,100) , ylim=c(0,3))
the range of values on the x-axis (“xlim=c(0,100)”) is chosen manually for this specific covariate; of course it is possible to use something like ylim = c( 0 , max(MyData$EF) ).

Now plot the fitted smoothing spline using the lines() function:
lines( sm.spline(MyData$EF , exp(predicted$fit)) , col = "red" , lwd = 0.8)
Note that the term prediction gives log-hazard-ratios; therefore exp(predicted$fit) is plotted against the values of the covariate. The sm.spline() function is necessary since the points of the plot appear in random order and density, according to the underlying dataset; a plain lines() function would produce just a chaotic pattern. Alternative:
plot(MyData$EF , exp(predicted$fit) , col = "red" , cex = 0.2)
produces a scattered plot that reflects the distribution of the underlying data – I do prefer adding a rug-plot on the bottom of the graph to illustrate this (see under).

… upper and lower confidence limits with dashed thinner lines

lines(sm.spline(MyData$EF , exp(predicted$fit + 1.96 * predicted$se)) , col = "orange" , lty = 2 , lwd = 0.4)
and
lines(sm.spline(MyData$EF , exp(predicted$fit - 1.96 * predicted$se)) , col = "orange" , lty = 2 , lwd = 0.4)

… a tiny horizontal line at hazard level 1, do see where the confidence limits cross:
abline( h = 1 , col = "lightgrey" , lty = 2 , lwd = 0.4)

… tiny tickmarks on the x-axes to reflect the distribution of the underlying data:
axis( side = 1 , at = MyData$EF, labels = F , tick = T , tcl = 0.4 , lwd.ticks = 0.1)

… and some fancy red tickmarks to mark minimum, lower hinge, median, upper hinge and maximum of the covariate in the dataset:
axis( side = 1 , at = fivenum(MyData$EF), labels = F , tick = T , tcl = -0.2 , lwd.ticks = 1 , col.ticks = "red")

Fancy customized smoothing spline fitted to the functional form of a covariate in a additive proportional hazard model
Fancy customized smoothing spline fitted to the functional form of a covariate in a additive proportional hazard model

Thats it!

4b) The easy way (works ONLY with MORE then 1 continous covariate) – predicting the terms can be omitted:

termplot(pham.fit, se=T, rug=T)

Resulting in …

The default termplot method for fitted smoothing splines
The default termplot method for fitted smoothing splines