Alternate R Code

This page contains R code that supplements or replaces code in the first printing of the second edition of RMS.

Examples of Replacement Code for Code That No Longer Works Because of Changes to R or R Packages

P. 516 Code to Produce Figure 20.7

## Simulate a proportional hazards situation with a quadratic covarate
## effect, and estimate the transformation with a linear fit.  Show
## smoothed martingale residuals to estimate true fit

require(rms)
n <- 1000
set.seed(1)
age <- 50 + 12 * rnorm(n)
cens <- 15 * runif(n)
h <- 0.02 * exp(0.04 * (age - 50) + 0.004 * (age - 50) ^ 2)
y <- -log(runif(n)) / h
e <- ifelse(y <= cens, 1, 0)
t <- pmin(y, cens)
S <- Surv(t, e)
f <- cph(S ~ age, x=TRUE, y=TRUE)
res <- resid(f)
ggplot(data.frame(age, res), aes(age, res)) + geom_smooth(method='loess') +
  ylab('Martingale Residual')

## Can also use ols to show the residual trend:
g <- ols(res ~ rcs(age, 5))
ggplot(Predict(g, age=10:100))
</highlight>

P. 516 Code to Produce Figure 20.10

## In this simulated example, a linear model is used and proportional
## hazards by definition is not satisfied

require(rms)
set.seed(1)
n <- 1000
x <- runif(n)
y <- runif(n) + 2 * x 
Srv <- Surv(y)
cox <- cph(Srv ~ x, x=TRUE, y=TRUE)
cox.zph(cox, "rank")   # Test for PH for each column of X
res <- resid(cox, "scaledsch")

ggplot(data.frame(y, res), aes(y, res)) + geom_smooth(method='loess') +
  geom_hline(yintercept = coef(cox)) +
    xlab('Time') + ylab('Scaled Schoenfeld Residual')