inst/doc/concordance.R

In survival: Survival Analysis

### R code from vignette source 'concordance.Rnw'

###################################################
### code chunk number 1: concordance.Rnw:23-32
###################################################
options(continue="  ", width=60)
options(SweaveHooks=list(fig=function() par(mar=c(4.1, 4.1, .3, 1.1))))
pdf.options(pointsize=10) #text in graph about the same as regular text
options(contrasts=c("contr.treatment", "contr.poly")) #ensure default

#require("survival")
#library(survival)
library(survival)
library(splines)


###################################################
### code chunk number 2: examples1a
###################################################
# x1 and y2 are both continuous variables
concordance(y2 ~ x1, data= anscombe)


###################################################
### code chunk number 3: examples1b
###################################################
# Fisher's iris data
fit1 <- glm(Species=="versicolor" ~ ., family=binomial, data=iris) 
concordance(fit1)  # equivalent to an AUC


###################################################
### code chunk number 4: examples1c
###################################################
# Anscombe data (all variables are continuous)
fit2 <- lm(y2 ~ x1 + x4, data= anscombe)
concordance(fit2)  # C
sqrt(summary(fit2)$r.squared) # R


###################################################
### code chunk number 5: examples1d
###################################################
# parametric survival 
fit3 <- survreg(Surv(time, status) ~ karno + age + trt, data=veteran)
concordance(fit3)


###################################################
### code chunk number 6: examples1e
###################################################
# 3 Cox models
fit4 <- coxph(Surv(time, status) ~ karno + age + trt, data=veteran)
fit5 <- update(fit4, . ~ . + celltype)
fit6 <- update(fit5, . ~ . + prior)
ctest <- concordance(fit4, fit5, fit6)
ctest


###################################################
### code chunk number 7: balance
###################################################
getOption("SweaveHooks")[["fig"]]()
# The balance figure for the concordance document
btree <- function(n) {
   tfun <- function(n, id, power) {
       if (n==1) id
       else if (n==2) c(2*id, id)
       else if (n==3) c(2*id, id, 2*id+1)
       else {
           nleft <- if (n== power*2) power  else min(power-1, n-power/2)
           c(tfun(nleft, 2*id, power/2), id,
             tfun(n-(nleft+1), 2*id +1, power/2))
           }
       }
   tfun(n, 1, 2^(floor(logb(n-1,2))))
   }

temp <- c(1,2,6,8, 9,12,14, 18, 19, 21, 23, 24, 27    )
indx <- btree(13)

xpos <- 1:15
xpos[4:7] <- tapply(xpos[8:15], rep(1:4, each=2), mean)
xpos[2:3] <- tapply(xpos[4:7], rep(1:2, each=2),mean)
xpos[1] <- mean(xpos[2:3])
ypos <- rep(4:1, c(1,2,4,8))

oldpar <- par(mar=c(1,1,1,1))

plot(xpos, ypos, type='n', xaxt='n', yaxt='n', bty='n',
     xlab="", ylab="")
temp2 <-  c(13,7,5,3,3,3,1,1,1,1,1,1,1)
#text(xpos[indx], ypos[indx], paste(temp, " (", temp2[indx], ")",  sep=''))
text(xpos[indx], ypos[indx], as.character(temp))

delta=.1
for (i in 1:6) {
    segments(xpos[i]-delta, ypos[i]-delta,
             xpos[2*i]+delta, ypos[2*i]+delta)
    segments(xpos[i]+delta, ypos[i]-delta,
             xpos[2*i+1]-delta, ypos[2*i+1] +delta)
}
par(oldpar)


###################################################
### code chunk number 8: concordance.Rnw:311-316
###################################################
# Concordance for a coxph object
concordance(fit4)

# Concordance using predictions from a Cox model
concordance(Surv(time, status) ~ predict(fit4), data = veteran, reverse = TRUE)


###################################################
### code chunk number 9: veteran2
###################################################
fit4b <- coxph(formula = Surv(time, status) ~ karno + age + trt + 
                   strata(celltype), data = veteran)
concordance(fit4b)
table(veteran$celltype)


###################################################
### code chunk number 10: amlexample
###################################################
afit <- survfit(Surv(time, status) ~1, aml, se = FALSE)
summary(afit, times = afit$time[1:6], censor = TRUE)


###################################################
### code chunk number 11: tmwt
###################################################
getOption("SweaveHooks")[["fig"]]()
colonfit <- coxph(Surv(time, status) ~ rx + nodes + extent, data = colon,
                 subset = (etype == 2))   # death only
cord1 <- concordance(colonfit, timewt="n",    ranks=TRUE)
cord2 <- concordance(colonfit, timewt="S",    ranks=TRUE)
cord3 <- concordance(colonfit, timewt="S/G",  ranks=TRUE)
cord4 <- concordance(colonfit, timewt="n/G2", ranks=TRUE)
temp <- c("n(t)"= coef(cord1), S=coef(cord2), "S/G"= coef(cord3), 
          "n/G2"= coef(cord4))
round(temp,5) # 4 different concordance estimates

# Plot the weights over time using the first 3 approaches
matplot(cord1$ranks$time/365.25, cbind(cord1$ranks$timewt,
                                       cord2$ranks$timewt,
                                       cord3$ranks$timewt),
        type= "l", lwd=2, col=c(1,2,4),
        xlab="Years since enrollment", ylab="Weight")
legend(1, 3000, c("n(t)", "nS(t-)", "nS(t-)/G(t-)"), lwd=2,
       col=c(1,2,4), lty=1:3, bty="n")

# Note that n/G2 and S/G are identical
all.equal(cord3$ranks$timewt,cord4$ranks$timewt)


###################################################
### code chunk number 12: manycurve
###################################################
getOption("SweaveHooks")[["fig"]]()
duo <- function(time, status, name, conf.int=FALSE) {
    sfit <- survfit(Surv(time, status) ~1)
    gfit <- survfit(Surv(time, max(status)-status) ~1)
    plot(sfit, conf.int=conf.int, xlab=name, lwd=2)
    lines(gfit, col=2, lwd=2, conf.int = conf.int)
}
oldpar <- par(mfrow=c(3,3), mar=c(5,5,1,1))
with(subset(colon, etype==1), duo(time/365.25, status, "NCCTG colon cancer"))
duo(flchain$futime/365.25, flchain$death, "Free light chain")
duo(kidney$time/12, kidney$status, "McGilchrist kidney")
duo(lung$time/365.25, lung$status, "advanced lung cancer")
duo(mgus2$futime/12, mgus2$death, "MGUS")
duo(nafld1$futime/365.25, nafld1$status, "NAFLD")

duo(pbc$time/365.25, pmin(pbc$status,1), "PBC")
with(rotterdam, duo(pmin(rtime, dtime)/365.25, pmax(recur, death), "Rotterdam"))

nfit <- coxph(Surv(futime/365.25, status) ~ age + male, data = nafld1)
znfit <- cox.zph(nfit, transform = 'identity')
plot(znfit[1], resid = FALSE)
par(oldpar)


###################################################
### code chunk number 13: nafld1
###################################################
nfit <- coxph(Surv(futime/365.25, status) ~ age + male, data = nafld1)
ncord1 <- concordance(nfit, timewt = "n")
ncord2 <- concordance(nfit, timewt = "S")
ncord3 <- concordance(nfit, timewt = "S/G")
ncord4 <- concordance(nfit, timewt = "n/G2")
temp <- c(n = coef(ncord1), S = coef(ncord2), 
          "S/G" = coef(ncord3), "n/G2" = coef(ncord4))
round(temp,6)


###################################################
### code chunk number 14: rankresid (eval = FALSE)
###################################################
## # code for Figure 3
## # pick a dataset with a smaller number of points, and non PH
## vfit <- coxph(Surv(time/365.25, status) ~ age + karno, data = veteran)
## temp <- concordance(vfit, ranks=TRUE)$rank
## # Two outliers at 999 days = 2.7 years stretch the axis too far
## plot(rank ~ time, data=temp, xlim=c(0,1.6),
##      xlab="Years", ylab="Rank residual")
## lines(lowess(temp$time, temp$rank, iter=1), lwd=2, col=2)
## abline(0, 0, lty=3)


###################################################
### code chunk number 15: rankresid2
###################################################
getOption("SweaveHooks")[["fig"]]()
# code for Figure 3
# pick a dataset with a smaller number of points, and non PH
vfit <- coxph(Surv(time/365.25, status) ~ age + karno, data = veteran)
temp <- concordance(vfit, ranks=TRUE)$rank
# Two outliers at 999 days = 2.7 years stretch the axis too far
plot(rank ~ time, data=temp, xlim=c(0,1.6),
     xlab="Years", ylab="Rank residual")
lines(lowess(temp$time, temp$rank, iter=1), lwd=2, col=2)
abline(0, 0, lty=3)


###################################################
### code chunk number 16: rotterdam (eval = FALSE)
###################################################
getOption("SweaveHooks")[["fig"]]()
## # should this be included?
## # recurrence free survival = earlier of recurrence and death
## rdata <- rotterdam
## rdata$rfs <-     with(rdata, ifelse(recur==1, 1, death))
## rdata$rfstime <- with(rdata, ifelse(recur==1, rtime, dtime))/ 365.25
## 
## rfit <- coxph(Surv(rfstime, rfs) ~ age + meno + grade + pspline(nodes), rdata)
## ctemp <- matrix(0, 100, 2)  # concordance and std err
## ctime <- seq(.1, 10, length=100)
## for (i in 1:100) { 
##     temp <- concordance(rfit, ymax=ctime[i])
##     ctemp[i,] <- c(temp$concordance, sqrt(temp$var))
## }
## yhat <- ctemp[,1] + outer(ctemp[,2], c(0, -1.96, 1.96), '*')
## matplot(ctime, yhat, type='l', lty=c(1,2,2), lwd=c(2,1,1), col=1, 
##         xlab="Upper cutoff", ylab="C", ylim=c(0.5,1))


###################################################
### code chunk number 17: test
###################################################
ctest <- concordance(fit4, fit5, fit6)
ctest

# compare concordance values of fit4 and fit5
contr <- c(-1, 1, 0)
dtest <- contr %*% coef(ctest)
dvar  <- contr %*% vcov(ctest) %*% contr

c(contrast=dtest, sd=sqrt(dvar), z=dtest/sqrt(dvar))


###################################################
### code chunk number 18: Cztrans
###################################################
zci <- function(fit, p=.95) {
    ilogist   <- function(p) log(p/(1-p))         # inverse logistic
    logistic  <- function(x) exp(x)/(1 + exp(x))  
    temp <- concordance(fit, influence =1)
    cminus <- temp$concordance - temp$dfbeta # values of concordance, without i

    newd    <- ilogist(temp$concordance) - ilogist(cminus) # dfbeta on new scale
    new.sd  <- sqrt(sum(newd^2))
    old.sd  <- sqrt(sum(temp$dfbeta^2))  # same as sqrt(temp$var)
    
    z <-  qnorm((1-p)/2)
    old.ci  <- temp$concordance + c(z, -z)*old.sd
    new.ci  <- logistic(ilogist(temp$concordance) + c(z, -z)* new.sd)
    rbind(old = old.ci, new= new.ci)
}

round(zci(colonfit), 4)


###################################################
### code chunk number 19: close
###################################################
set.seed(1953)
ytest <- matrix(rexp(20), ncol=2) %*% chol(matrix(c(1, .98, .98, 1), 2))
cor(ytest)
lfit <- lm(ytest[,1] ~ ytest[,2])
zci(lfit)