Assessment of Discrimination in Survival Analysis (C-statistics, etc)

References

Discrimination

Different C-statistics and related measures implemented in R

See individual examples below for links to the original papers.

For the explanation of the difference between cumulative/dynamic AUCs vs incident/dynamic AUCs, see the paper by Heargery et al ( http://www.statmed.medicina.unimib.it/statisticalps2011/materiale/Heagerty%20and%20Zheng,%20Biometrics%202005.pdf ).

As I understand it,

Thus, the idea behind incident/dynamic AUCs is closer to the idea of hazard (dynamically changing instantaneous incidence at a given time), and it can handle time-varying predictors. But the question answered by cumulative/dynamic AUCs may be more clinically relevant, e.g., does this model discriminate if I will survive next five years .

CRAN Task View: Survival Analysis: Predictions and Prediction Performance

PBC dataset

       age:       in years
       albumin:   serum albumin (g/dl)
       alk.phos:  alkaline phosphotase (U/liter)
       ascites:   presence of ascites
       ast:       aspartate aminotransferase, once called SGOT (U/ml)
       bili:      serum bilirunbin (mg/dl)
       chol:      serum cholesterol (mg/dl)
       copper:    urine copper (ug/day)
       edema:     0 no edema, 0.5 untreated or successfully treated
                  1 edema despite diuretic therapy
       hepato:    presence of hepatomegaly or enlarged liver
       id:        case number
       platelet:  platelet count
       protime:   standardised blood clotting time
       sex:       m/f
       spiders:   blood vessel malformations in the skin
       stage:     histologic stage of disease (needs biopsy)
       status:    status at endpoint, 0/1/2 for censored, transplant, dead
       time:      number of days between registration and the earlier of death,
                  transplantion, or study analysis in July, 1986
       trt:       1/2/NA for D-penicillmain, placebo, not randomised
       trig:      triglycerides (mg/dl)

Load the PBC dataset and modify for later use

library(survival)
data(pbc)

pbc <- within(pbc, {
    ## transplant (1) and death (2) are considered events, and marked 1
    event <- as.numeric(status %in% c(1,2))

    ## Create a survival vector
    Surv <- Surv(time, event)
})

Kaplan-Meier plots

par(mar = c(2,2,2,2))
layout(matrix(1:4, byrow = T, ncol = 2))
library(rms)
survplot(survfit(Surv ~ 1, data = pbc))
survplot(survfit(Surv ~ age >= median(age), data = pbc),         label.curves = list(method = "arrow", cex = 1.2))
survplot(survfit(Surv ~ sex, data = pbc),                        label.curves = list(method = "arrow", cex = 1.2))
survplot(survfit(Surv ~ albumin >= median(albumin), data = pbc), label.curves = list(method = "arrow", cex = 1.2))

plot of chunk unnamed-chunk-3

AUC by logistic regression models

Completely ignore the time variable

It is the simplest method. Logistic regression is used instead of Cox regression model.

Completely ignore the time variable and use the outcome variable as a binary outcome variable. It does not take into acount the variable length of follow-up.

## Load epicalc package to calcuate AUC
library(epicalc)

## Model with age and sex
logit.age.sex <- glm(event ~ age + sex, data = pbc, family = binomial)
lroc(logit.age.sex, graph = F)$auc
[1] 0.5879
## Model with age, sex, and albumin
logit.age.sex.albumin <- glm(event ~ age + sex + albumin, data = pbc, family = binomial)
lroc(logit.age.sex.albumin, graph = F)$auc
[1] 0.6751

Cut the follow up at a specifict time point

Only events that occured within two years are considered events and others are treated as non-events. This method can be valid if the specified time is short enough so that there are few censored subjects. Most people have complete follow-up in this situation.

## Create a variable indicating 2-year event
pbc <- within(pbc, {
    outcome2yr <- NA
    outcome2yr[(event == 1) & (time <= 2 * 365)] <- 1 # event+ within two years
    outcome2yr[(event == 0) | (time  > 2 * 365)] <- 0 # otherwise
})

## 2-year outcome model with age and sex
logit.age.sex <- glm(outcome2yr ~ age + sex, data = pbc, family = binomial)
lroc(logit.age.sex, graph = F)$auc
[1] 0.6484
## 2-year outcome model with age, sex, and albumin
logit.age.sex.albumin <- glm(outcome2yr ~ age + sex + albumin, data = pbc, family = binomial)
lroc(logit.age.sex.albumin, graph = F)$auc
[1] 0.7918

Fit Cox regression models for later use

The linear predictors (lp), when exponetiated, will provide the predicted hazard ratios for individuals. Thus these can be used as the summary predictors calculated from multiple raw predictors.

## Null model
coxph.null             <- coxph(Surv ~ 1, data = pbc)
coxph.null
Call:  coxph(formula = Surv ~ 1, data = pbc)

Null model
  log likelihood= -1009 
  n= 418 
## Model with age and sex
coxph.age.sex          <- coxph(Surv ~ age + sex, data = pbc)
coxph.age.sex
Call:
coxph(formula = Surv ~ age + sex, data = pbc)

        coef exp(coef) se(coef)     z      p
age   0.0221     1.022  0.00727  3.04 0.0024
sexf -0.3000     0.741  0.20975 -1.43 0.1500

Likelihood ratio test=12  on 2 df, p=0.0025  n= 418, number of events= 186 
## Model with age, sex, and albumin
coxph.age.sex.albumin  <- coxph(Surv ~ age + sex + albumin, data = pbc)
coxph.age.sex.albumin
Call:
coxph(formula = Surv ~ age + sex + albumin, data = pbc)

           coef exp(coef) se(coef)     z       p
age      0.0148     1.015  0.00749  1.97 4.9e-02
sexf    -0.3467     0.707  0.21079 -1.64 1.0e-01
albumin -1.4066     0.245  0.17096 -8.23 2.2e-16

Likelihood ratio test=73.9  on 3 df, p=6.66e-16  n= 418, number of events= 186 
## These models are significantly different by likelihood ratio test
anova(coxph.age.sex, coxph.age.sex.albumin, test = "LRT")
Analysis of Deviance Table
 Cox model: response is  Surv
 Model 1: ~ age + sex
 Model 2: ~ age + sex + albumin
  loglik Chisq Df P(>|Chi|)    
1  -1003                       
2   -972  61.9  1   3.7e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## Put linear predictors ("lp") into pbc dataset
pbc$lp.null            <- predict(coxph.null, type = "lp")
pbc$lp.age.sex         <- predict(coxph.age.sex, type = "lp")
pbc$lp.age.sex.albumin <- predict(coxph.age.sex.albumin, type = "lp")

Harrell’s C or concordance

Trying to compare (test) the difference in C is not recommended.

http://stats.stackexchange.com/questions/17480/how-to-do-roc-analysis-in-r-with-a-cox-model/17517#17517

Hmisc::rcorrcens

A larger marker value is considered to be associated with a longer survival by this function. Thus, the linear predictor (the higher, the worse) needs to be negated

library(Hmisc)

## Model with age and sex
rcorrcens(formula = Surv ~ I(-1 * lp.age.sex), data = pbc)

Somers' Rank Correlation for Censored Data    Response variable:Surv

                       C   Dxy  aDxy    SD    Z      P   n
I(-1 * lp.age.sex) 0.569 0.138 0.138 0.043 3.18 0.0015 418
## Model with age, sex, and albumin
rcorrcens(formula = Surv ~ I(-1 * lp.age.sex.albumin), data = pbc)

Somers' Rank Correlation for Censored Data    Response variable:Surv

                               C   Dxy  aDxy    SD    Z P   n
I(-1 * lp.age.sex.albumin) 0.703 0.407 0.407 0.041 9.86 0 418

survival package

Actually, the summary method for coxph objects prints “Concordance” (five lines from bottom), which is the same thing as the Harrells’C, and \( R^2 \).

## Model with age and sex
summary(coxph.age.sex)
Call:
coxph(formula = Surv ~ age + sex, data = pbc)

  n= 418, number of events= 186 

         coef exp(coef) se(coef)     z Pr(>|z|)   
age   0.02210   1.02234  0.00727  3.04   0.0024 **
sexf -0.29995   0.74085  0.20975 -1.43   0.1527   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

     exp(coef) exp(-coef) lower .95 upper .95
age      1.022      0.978     1.008      1.04
sexf     0.741      1.350     0.491      1.12

Concordance= 0.569  (se = 0.023 )
Rsquare= 0.028   (max possible= 0.992 )
Likelihood ratio test= 12  on 2 df,   p=0.0025
Wald test            = 12.2  on 2 df,   p=0.0022
Score (logrank) test = 12.3  on 2 df,   p=0.00214
## Model with age, sex, and albumin
summary(coxph.age.sex.albumin)
Call:
coxph(formula = Surv ~ age + sex + albumin, data = pbc)

  n= 418, number of events= 186 

            coef exp(coef) se(coef)     z Pr(>|z|)    
age      0.01476   1.01487  0.00749  1.97    0.049 *  
sexf    -0.34672   0.70700  0.21079 -1.64    0.100 .  
albumin -1.40664   0.24496  0.17096 -8.23  2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

        exp(coef) exp(-coef) lower .95 upper .95
age         1.015      0.985     1.000     1.030
sexf        0.707      1.414     0.468     1.069
albumin     0.245      4.082     0.175     0.342

Concordance= 0.703  (se = 0.023 )
Rsquare= 0.162   (max possible= 0.992 )
Likelihood ratio test= 73.9  on 3 df,   p=6.66e-16
Wald test            = 83.1  on 3 df,   p=0
Score (logrank) test = 82.2  on 3 df,   p=0

Gonen and Heller Concordance Index for Cox models

CPE package

library(CPE)
## Model with age and sex
phcpe(coxph.age.sex, CPE.SE = TRUE)
$CPE
[1] 0.5724

$CPE.SE
[1] 0.02023
## Model with age, sex, and albumin
phcpe(coxph.age.sex.albumin, CPE.SE = TRUE)
$CPE
[1] 0.6615

$CPE.SE
[1] 0.01531

clinfun package

library(clinfun)
## Model with age and sex
coxphCPE(coxph.age.sex)
       CPE smooth.CPE     se.CPE 
   0.57244    0.57236    0.02023 
## Model with age, sex, and albumin
coxphCPE(coxph.age.sex.albumin)
       CPE smooth.CPE     se.CPE 
   0.66150    0.66128    0.01531 

Cumulative case/dynamic control ROC/AUC using survivalROC package (Heagerty, 2000)

library(survivalROC)

## Define a function
fun.survivalROC <- function(lp, t) {
    res <- with(pbc,
                survivalROC(Stime        = time,
                            status       = event,
                            marker       = get(lp),
                            predict.time = t,
                            method       = "KM"))       # KM method without smoothing

    ## Plot ROCs
    with(res, plot(TP ~ FP, type = "l", main = sprintf("t = %.0f, AUC = %.2f", t, AUC)))
    abline(a = 0, b = 1, lty = 2)

    res
}

## 2 x 5 layout
layout(matrix(1:10, byrow = T, ncol = 5))

## Model with age and sex
res.survivalROC.age.sex <- lapply(1:10 * 365.25, function(t) {
    fun.survivalROC(lp = "lp.age.sex", t)
})

plot of chunk unnamed-chunk-11


## Model with age, sex, and albumin
res.survivalROC.age.sex.albumin <- lapply(1:10 * 365.25, function(t) {
    fun.survivalROC(lp = "lp.age.sex.albumin", t)
})

plot of chunk unnamed-chunk-11

Incident case/dynamic control ROC/AUC using risksetROC package (Heagerty, 2005)

1- to 10-year ROCs

library(risksetROC)

## Define a function
fun.risksetROC <- function(lp, t) {
    res <- with(pbc,
                risksetROC(Stime        = time,
                           status       = event,
                           marker       = get(lp),
                           predict.time = t,
                           plot         = FALSE))

    ## Plot ROCs
    with(res, plot(TP ~ FP, type = "l", main = sprintf("t = %.0f, AUC = %.2f", t, AUC)))
    abline(a = 0, b = 1, lty = 2)

    res
}

## 2 x 5 layout
layout(matrix(1:10, byrow = T, ncol = 5))

## Model with age and sex
risksetROC.age.sex <- lapply(365.25 * 1:10, function(t) {
    fun.risksetROC(lp = "lp.age.sex", t)
})

plot of chunk unnamed-chunk-12


## Model with age, sex, and albumin
risksetROC.age.sex.albumin <- lapply(365.25 * 1:10, function(t) {
    fun.risksetROC(lp = "lp.age.sex.albumin", t)
})

plot of chunk unnamed-chunk-12

Up to 10-year AUCs

## 1 x 2 layout
layout(matrix(1:2, byrow = T, ncol = 2))

## Model with age and sex
risksetAUC.age.sex <- with(pbc,
                           risksetAUC(Stime        = time,
                                      status       = event,
                                      marker       = lp.age.sex,
                                      tmax         = 10 * 365.25))
title(sprintf("t = %.0f, iAUC = %.2f", 10 * 365.25, risksetAUC.age.sex$Cindex))

## Model with age, sex, and albumin
risksetAUC.age.sex.albumin <- with(pbc,
                           risksetAUC(Stime        = time,
                                      status       = event,
                                      marker       = lp.age.sex.albumin,
                                      tmax         = 10 * 365.25))
title(sprintf("t = %.0f, iAUC = %.2f", 10 * 365.25, risksetAUC.age.sex.albumin$Cindex))

plot of chunk unnamed-chunk-13

Various methods provided survAUC package (Potapov et al)

These can calculate multiple time-dependent ROC at once, and also compute summary measures of a time-dependent AUC curve (iAUC)

These need a training dataset and a test dataset. The same data can be given to both, and it works although I am not sure if this is correct.

Functions

Time-dependent AUCs for the age sex model are calculated by various methods. (1- to 10-year AUCs)

library(survAUC)

## *Cumulative case*/dynamic control AUC by Chambless and Diao (Stat Med 2006;25:3474-3486.)
res.AUC.cd <- AUC.cd(Surv.rsp     = pbc$Surv,
                     Surv.rsp.new = pbc$Surv,
                     lp           = coxph.age.sex$linear.predictors,
                     lpnew        = coxph.age.sex$linear.predictors,
                     times        = 1:10 * 365.25)
res.AUC.cd$iauc
[1] 0.5895
plot(res.AUC.cd)

plot of chunk unnamed-chunk-14


## *Cumulative case*/dynamic control AUC by Hung and Chiang (Can J Stat 2010;38:8-26)
res.AUC.hc <- AUC.hc(Surv.rsp     = pbc$Surv,
                     Surv.rsp.new = pbc$Surv,
                     ## lp           = coxph.age.sex$linear.predictors,
                     lpnew        = coxph.age.sex$linear.predictors,
                     times        = 1:10 * 365.25)
res.AUC.hc$iauc
[1] 0.6005
plot(res.AUC.hc)

plot of chunk unnamed-chunk-14


## *Incident case* or *Cumulative case*/dynamic control AUC by Song and Zhou (Biometrics 2011;67:906-16)
res.AUC.sh <- AUC.sh(Surv.rsp     = pbc$Surv,
                     Surv.rsp.new = pbc$Surv,
                     lp           = coxph.age.sex$linear.predictors,
                     lpnew        = coxph.age.sex$linear.predictors,
                     times        = 1:10 * 365.25)
res.AUC.sh$iauc
[1] 0.573
plot(res.AUC.sh)

plot of chunk unnamed-chunk-14


## *Cumulative case*/dynamic control AUC by Uno et al.
## (http://biostats.bepress.com/cgi/viewcontent.cgi?article=1041&context=harvardbiostat)
res.AUC.uno <- AUC.uno(Surv.rsp     = pbc$Surv,
                       Surv.rsp.new = pbc$Surv,
                       ## lp           = coxph.age.sex$linear.predictors,
                       lpnew        = coxph.age.sex$linear.predictors,
                       times        = 1:10 * 365.25)
res.AUC.uno$iauc
[1] 0.5811
plot(res.AUC.uno)

plot of chunk unnamed-chunk-14

Summary measures (10 years when applicable)

## C-statistic by Uno et al. (Stat Med 2011;30:1105-1117)
res.UnoC <- UnoC(Surv.rsp     = pbc$Surv,
                 Surv.rsp.new = pbc$Surv,
                 ## lp           = coxph.age.sex$linear.predictors,
                 lpnew        = coxph.age.sex$linear.predictors,
                 time        = 10 * 365.25)       # Upper limit, not time points
res.UnoC
[1] 0.5608

## C-statistic by Begg et al. (Stat Med 2000;19:1997-2014)
res.BeggC <- BeggC(Surv.rsp     = pbc$Surv,
                   Surv.rsp.new = pbc$Surv,
                   lp           = coxph.age.sex$linear.predictors,
                   lpnew        = coxph.age.sex$linear.predictors)
res.BeggC
[1] 0.5519

## Gonen and Heller’s Concordance Index for Cox PH models (Biometrika 2005;92:965-970)
res.GHCI <- GHCI(## Surv.rsp     = pbc$Surv,
                 ## Surv.rsp.new = pbc$Surv,
                 ## lp           = coxph.age.sex$linear.predictors,
                 lpnew        = coxph.age.sex$linear.predictors)
res.GHCI
[1] 0.5717

*R2-type Coefficients for Cox models *

## O'Quigley et al. (Stat Med 2005;24:479-489)
res.OXS <- OXS(Surv.rsp = pbc$Surv,
               lp       = coxph.age.sex$linear.predictors,
               lp0      = coxph.null$linear.predictors)
res.OXS
[1] 0.06092

## Nagelkerke (Biometrika 1991;78:691-692)
res.Nagelk <- Nagelk(Surv.rsp = pbc$Surv,
                     lp       = coxph.age.sex$linear.predictors,
                     lp0      = coxph.null$linear.predictors)
res.Nagelk
[1] 0.02779

## Xu et al. (Journal of Nonparametric Statistics 1999;12:83-107)
res.XO <- XO(Surv.rsp = pbc$Surv,
             lp       = coxph.age.sex$linear.predictors,
             lp0      = coxph.null$linear.predictors)
res.XO
[1] 0.04166

Uno methods for C-statistics and IDI/NRI implemented in survC1 and survIDINRI

## Create numeric variable for sex
pbc$female <- as.numeric(pbc$sex == "f")

C-statistics (10-years follow up) using survC1 package

library(survC1)
## C-statistic for age sex model
unoC.age.female         <- Est.Cval(mydata = pbc[,c("time","event","age","female")],           tau = 10 * 365.25)
unoC.age.female$Dhat
[1] 0.5618
## C-statistic for age sex albumin model
unoC.age.female.albumin <- Est.Cval(mydata = pbc[,c("time","event","age","female","albumin")], tau = 10 * 365.25)
unoC.age.female.albumin$Dhat
[1] 0.672

Comparison of C-statistics

  mydata: Input data. The 1st column should be time-to-event, and the
          2nd column is event indicator (1=event, 0=censor).
   covs0: A matrix that consists of a set of predictors for a base
          model (Model 0)
   covs1: A matrix that consists of a set of predictors for a new model
          (Model 1)
     tau: Truncation time. The resulting C tells how well the given
          prediction model works in predicting events that occur in the
          time range from 0 to ‘tau’. Note that the survival function
          for the underlying censoring time distribution needs to be
          positive at ‘tau’.
     itr: Iteration of perturbation-resampling.
## Comaprison of C-statistics
uno.C.delta <- Inf.Cval.Delta(mydata = pbc[,c("time","event")],
                              covs0  = pbc[,c("age","female")],           # age sex model
                              covs1  = pbc[,c("age","female","albumin")], # age sex albumin model
                              tau    = 10 * 365.25,                       # Trucation time (max time to consider)
                              itr = 1000)                                 # Iteration of perturbation-resampling
uno.C.delta
          Est      SE Lower95 Upper95
Model1 0.6720 0.02137 0.63009  0.7138
Model0 0.5618 0.02266 0.51736  0.6061
Delta  0.1102 0.02367 0.06384  0.1566

IDI, continous NRI, and median improvement (10-years follow up) using survIDINRI

library(survIDINRI)
res.IDI.INF <- IDI.INF(indata = pbc[,c("time","event")],
                       covs0 = pbc[,c("age","female")],           # age sex model
                       covs1 = pbc[,c("age","female","albumin")], # age sex albumin model
                       t0 = 10 * 365.25,
                       npert = 300, npert.rand = NULL, seed1 = NULL, alpha = 0.05)

## M1 IDI; M2 continuous NRI; M3 median improvement
IDI.INF.OUT(res.IDI.INF)
    Est. Lower Upper
M1 0.101 0.027 0.162
M2 0.254 0.015 0.453
M3 0.094 0.018 0.180
## M1 red area; M2 distance between black points; M3 distance between gray points
IDI.INF.GRAPH(res.IDI.INF)

plot of chunk unnamed-chunk-20