Rforge

_

Plot Function For Additive Cox Regression

A Cox-regression model with R is achieved 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)

When I started out writing a plot function for the spline fit, the termplot() function failed for univariate models this is fixed now. So the main incentive to use plotHR() is gone.

On the other side there is a lot of customization possible with plotHR() which is not with termplot() and one my find the syntax interesting to see how the internals of the coxph() function work. I learned at least a lot from writing this.

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

Download plotHR() to your harddisk, remember the path and execute
source("C:Path/to/plotHR_0.7.R")
before using the function. Note the forward slash “/” in the path description in R – using the Windows backslash “\” will fail.

You might also find it useful to open the downloaded file plotHR_0.x.R and have a look inside. I am trying to keep the code as readable as possible and to comment extensively. Feel free to send feedback, bug reports, improvements and of course to addapt it to your needs.

If you need to remember the possible options in the function call after you have “sourced” the function you can also type
plotHR
in the R-command line without brackets and R will display the syntax of the function.

Download plotHR()

plotHR_0.9.R

Usage

plotHR( model , terms = 1 , se = TRUE , rug = "ticks" , 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.
• dens.pos = 1
• 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

Version history

Note: I have rewritten the function several times since I wrote the initial post.

• V0.9: Mainly rewrite of the density plot. rug = “density” will produce acceptable densityplots on the x-axis with marks for 2.5-, 25-, 50-, 75- and 97.5-percentiles. Also the confidence shade is marked at the same percentiles. This will only appear when rug = “density” is specified. The position of the densityplot in the main plot panel is still a concern: It should work for one-panel plots. In multipanel-plots the densityplot will move under the x-axis. This can be manually adjusted by selecting a fitting value for dens.pos to move it up again. For a 1×3-panel I had to use dens.pos = 0.6.
• V0.7 is a complete rewrite.
  • The local datasset is trimmed to include only model covariates and complete observations only. The uncomplete observations are omitted in the cox model anyway. This removes difficulties which arose from missing covariate values in the dataset.
  • x.log option is depretiated – it is planned to use the xlog option from plot.default() on the long run.
  • the rug option is extended to FALSE, “density” or “ticks” (the default). The option rug = “density” is still a bit alpha’ish but includes a density plot of the covariate distribution on the x-axis.
• 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).

Reference

The additive model and the plotting method is based on:

Therneau TM, Grambsch PM: Modeling Survival Data: Extending the Cox Model. Springer 2000.

The relevant part is on page 108ff, chapter 4, Modeling the functional form of covariates.

Therneau is the maintainer of the “survival” package in R. I wonder why the additive modelling of a covariate is not mentioned more explicitly in “library(survival)”.



Share this: