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