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))

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')
Topic revision: r2 - 23 Sep 2015, FrankHarrell
 

This site is powered by FoswikiCopyright © 2013-2022 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding Vanderbilt Biostatistics Wiki? Send feedback