How to plot a Cox hazard model with splines

Question

I have a following model:

coxph(Surv(fulength, mortality == 1) ~ pspline(predictor))

where is fulength is a duration of follow-up (including mortality), predictor is a predictor of mortality.

The output of the command above is this:

                         coef  se(coef) se2    Chisq DF   p    
pspline(predictor), line 0.174 0.0563   0.0562 9.52  1.00 0.002
pspline(predictor), nonl                       4.74  3.09 0.200

How can I plot this model so that I get the nice curvy line with 95% confidence bands and hazard ratio on the y axis? What I am aiming for is something similar to this:

Kiger

stat.ethz.ch/R-manual/R-devel/library/graphics/html/curve.html r-bloggers.com/plotting-95-confidence-bands-in-r-2 – Sugden

I use Frank Harrell's rms/Hmisc packages which are probably capable of delivering something very much like that output although I don't know about the right-hand plot. I'm statistically offended by plotting a normal distribution under the Male and Female results. I don't know if rms supports psplines because Frank prefers restricted cubic splines, but if you post some data I'd be happy to give it a try. – Antiphon

Thanks BondedDust. how did you install these packages? When I try to install, I get error messages like this installation of package ‘TH.data’ had non-zero exit status – Kiger

Assuming you have an open Internet connection and current version of R, then this should succeed: install.packages(c("Hmisc","rms"), dependencies=TRUE ). I do not recognize the package that is being reported as missing. I suspect a syntax error on your part. – Antiphon

Answer 1

This is when you get when you run the first example in ?cph of the rms-package:

n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n, 
              rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
dt <- -log(runif(n))/h
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
dd <- datadist(age, sex)
options(datadist='dd')
S <- Surv(dt,e)

f <- cph(S ~ rcs(age,4) + sex, x=TRUE, y=TRUE)
cox.zph(f, "rank")             # tests of PH
anova(f)
plot(Predict(f, age, sex)) # plot age effect, 2 curves for 2 sexes

Because the rms/Hmisc package combo uses lattice plots, annotation with a marginal age-density feature would need to be done with lattice-functions. On the other hand, if you want to change the response value to relative hazard you can just add a 'fun=exp' argument to the Predict call and relable the graph to get:

png(); plot(Predict(f, age, sex, fun=exp), ylab="Relative Hazard");dev.off()

Antiphon

I was able to install Hmisc but not rms. the error I get is this: ERROR: installing Rd objects failed for package ‘quantreg’ ERROR: dependency ‘TH.data’ is not available for package ‘multcomp’ ERROR: dependencies ‘quantreg’, ‘multcomp’ are not available for package ‘rms’ – Kiger

Sounds like the repository you are using has a broken or missing versions of multcomp and quantreg. I just tried installing the binary Mac version of multcomp from the Berkeley repo and had no errors. – Antiphon

Figured out the alternative way to install the missing packages. I get through all the steps but when I try to plot plot(Predict(f, predictor, fun=exp), ylab="Relative Hazard") I get the error Error in value.chk(at, which(name == n), NA, np, lim): variable predictor does not have limits defined by datadist – Kiger

Sounds like an error I often make: forgetting to set the datadist properly. You do need the lines dd <- datadist(age, sex); options(datadist='dd') or equivalent run on your data source. – Antiphon

I see, thanks for the tip. I fixed it. In the code you provided, cox.zph(f, "rank") # tests of PH what does this line mean? The way I understand this, p values for each slopes plus overall p value? – Kiger

cox.zph from the survival package returns tests for the proportionality assumption. Large p-values are desirable rather than small ones because one would hope that evidence against the assumption is weak. – Antiphon

I know this is an old post, but in the two curves above (male vs female), is there a way to obtain a p value of whether the two curves are different from each other?? – Kiger

It's really quite easy to obtain a p-value. Take the full model and compare to a reduced model where the sex term has been omitted. The difference in deviance ( -2*LL_full --2*LL_reduced) ) would be referenced against a chi-square value with one degree of freedom. If you cannot find a prior SO question (or a worked example using Google) showing how to do that, you should post another question. – Antiphon

Another follow-up, I haven't been able to find a good source of when to use the default "km" vs. "rank" vs. "identity", could you refer a link or if brief enough, describe here? – Scopophilia

Seems like another question, needing a minimal reproducible example. – Antiphon