# fit <- cph(Surv(start, stop, event) ~ (age + surgery)* transplant,
# data=jasa1)
# mresid <- resid(fit, collapse=jasa1$id)
# Get unadjusted relationships for several variables
# Pick one variable that's not missing too much, for fit
n <- 1000 # define sample size
set.seed(17) # so can reproduce the results
age <- rnorm(n, 50, 10)
blood.pressure <- rnorm(n, 120, 15)
cholesterol <- rnorm(n, 200, 25)
sex <- factor(sample(c('female','male'), n,TRUE))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
d.time <- -log(runif(n))/h
death <- ifelse(d.time <= cens,1,0)
d.time <- pmin(d.time, cens)
f <- cph(Surv(d.time, death) ~ age + blood.pressure + cholesterol, iter.max=0)
res <- resid(f) # This re-inserts rows for NAs, unlike f$resid
yl <- quantile(res, c(10/length(res),1-10/length(res)), na.rm=TRUE)
# Scale all plots from 10th smallest to 10th largest residual
par(mfrow=c(2,2), oma=c(3,0,3,0))
p <- function(x) {
s <- !is.na(x+res)
plot(lowess(x[s], res[s], iter=0), xlab=label(x), ylab="Residual",
ylim=yl, type="l")
}
p(age); p(blood.pressure); p(cholesterol)
mtext("Smoothed Martingale Residuals", outer=TRUE)
# Assess PH by estimating log relative hazard over time
f <- cph(Surv(d.time,death) ~ age + sex + blood.pressure, x=TRUE, y=TRUE)
r <- resid(f, "scaledsch")
tt <- as.numeric(dimnames(r)[[1]])
par(mfrow=c(3,2))
for(i in 1:3) {
g <- areg.boot(I(r[,i]) ~ tt, B=20)
plot(g, boot=FALSE) # shows bootstrap CIs
} # Focus on 3 graphs on right
# Easier approach:
plot(cox.zph(f)) # invokes plot.cox.zph
par(mfrow=c(1,1))
Run the code above in your browser using DataLab